Reference documentation for deal.II version GIT relicensing-539-geafd9ce159 2024-05-06 18:40: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\}}\)
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 active_fe_index = numbers::invalid_unsigned_int;
402
403 Assert(fe_collection->size() > 0,
404 ExcMessage("Incoming hp::FECollection can not be empty."));
405 Assert(
406 mapping_collection->size() == fe_collection->size() ||
407 mapping_collection->size() == 1,
409 "Size of hp::MappingCollection must be the same as hp::FECollection or 1."));
410 Assert(
411 q_collection.size() == fe_collection->size() || q_collection.size() == 1,
413 "Size of hp::QCollection<dim> must be the same as hp::FECollection or 1."));
414 Assert(
415 q_collection_1D.size() == fe_collection->size() ||
416 q_collection_1D.size() == 1,
418 "Size of hp::QCollection<1> must be the same as hp::FECollection or 1."));
419
420 // For each element in fe_collection, create ::FEInterfaceValues
421 // objects to use on the non-intersected cells.
422 fe_values_inside_full_quadrature.resize(fe_collection->size());
423 fe_values_outside_full_quadrature.resize(fe_collection->size());
424 for (unsigned int fe_index = 0; fe_index < fe_collection->size();
425 ++fe_index)
426 {
427 const unsigned int mapping_index =
428 mapping_collection->size() > 1 ? fe_index : 0;
429 const unsigned int q_index = q_collection.size() > 1 ? fe_index : 0;
430
431 fe_values_inside_full_quadrature[fe_index].emplace(
432 (*mapping_collection)[mapping_index],
433 (*fe_collection)[fe_index],
434 q_collection[q_index],
435 region_update_flags.inside);
436 fe_values_outside_full_quadrature[fe_index].emplace(
437 (*mapping_collection)[mapping_index],
438 (*fe_collection)[fe_index],
439 q_collection[q_index],
440 region_update_flags.outside);
441 }
442 }
443
444
445
446 template <int dim>
447 template <bool level_dof_access>
448 void
451 const unsigned int face_no,
452 const std::function<void(::FEInterfaceValues<dim> &)> &call_reinit)
453 {
454 current_face_location =
455 mesh_classifier->location_to_level_set(cell, face_no);
456 active_fe_index = cell->active_fe_index();
457
458 // These objects were created with a quadrature based on the previous cell
459 // and are thus no longer valid.
460 fe_values_inside.reset();
461 fe_values_outside.reset();
462
463 switch (current_face_location)
464 {
466 {
467 call_reinit(*fe_values_inside_full_quadrature.at(active_fe_index));
468 break;
469 }
471 {
472 call_reinit(*fe_values_outside_full_quadrature.at(active_fe_index));
473 break;
474 }
476 {
477 const unsigned int mapping_index =
478 mapping_collection->size() > 1 ? active_fe_index : 0;
479 const unsigned int q1D_index =
480 q_collection_1D.size() > 1 ? active_fe_index : 0;
481
482 face_quadrature_generator.set_1D_quadrature(q1D_index);
483 face_quadrature_generator.generate(cell, face_no);
484
485 const Quadrature<dim - 1> &inside_quadrature =
486 face_quadrature_generator.get_inside_quadrature();
487 const Quadrature<dim - 1> &outside_quadrature =
488 face_quadrature_generator.get_outside_quadrature();
489
490 // Even if a cell is formally intersected the number of created
491 // quadrature points can be 0. Avoid creating an FEInterfaceValues
492 // object if that is the case.
493 if (inside_quadrature.size() > 0)
494 {
495 fe_values_inside.emplace((*mapping_collection)[mapping_index],
496 (*fe_collection)[active_fe_index],
497 inside_quadrature,
498 region_update_flags.inside);
499
500 call_reinit(*fe_values_inside);
501 }
502
503 if (outside_quadrature.size() > 0)
504 {
505 fe_values_outside.emplace((*mapping_collection)[mapping_index],
506 (*fe_collection)[active_fe_index],
507 outside_quadrature,
508 region_update_flags.outside);
509
510 call_reinit(*fe_values_outside);
511 }
512 break;
513 }
514 default:
515 {
517 break;
518 }
519 }
520 }
521
522
523
524 template <int dim>
525 const std::optional<::FEInterfaceValues<dim>> &
527 {
528 if (current_face_location == LocationToLevelSet::inside)
529 return fe_values_inside_full_quadrature.at(active_fe_index);
530 else
531 return fe_values_inside;
532 }
533
534
535
536 template <int dim>
537 const std::optional<::FEInterfaceValues<dim>> &
539 {
540 if (current_face_location == LocationToLevelSet::outside)
541 return fe_values_outside_full_quadrature.at(active_fe_index);
542 else
543 return fe_values_outside;
544 }
545
546
547#include "fe_values.inst"
548
549} // namespace NonMatching
const std::optional<::FEInterfaceValues< dim > > & get_outside_fe_values() const
Definition fe_values.cc:538
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:655
typename FaceQuadratureGenerator< dim >::AdditionalData AdditionalData
Definition fe_values.h:492
const std::optional<::FEInterfaceValues< dim > > & get_inside_fe_values() const
Definition fe_values.cc:526
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
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:449
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:264
void push_back(const Quadrature< dim_in > &new_quadrature)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
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