Reference documentation for deal.II version GIT relicensing-249-g48dc7357c7 2024-03-29 12:30: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 {
153 current_cell_location = mesh_classifier->location_to_level_set(cell);
154 active_fe_index = cell->active_fe_index();
155
156 // These objects were created with a quadrature based on the previous cell
157 // and are thus no longer valid.
158 fe_values_inside.reset();
159 fe_values_surface.reset();
160 fe_values_outside.reset();
161
162 switch (current_cell_location)
163 {
165 {
166 fe_values_inside_full_quadrature.at(active_fe_index)->reinit(cell);
167 break;
168 }
170 {
171 fe_values_outside_full_quadrature.at(active_fe_index)->reinit(cell);
172 break;
173 }
175 {
176 const unsigned int mapping_index =
177 mapping_collection->size() > 1 ? active_fe_index : 0;
178
179 const unsigned int q1D_index =
180 q_collection_1D.size() > 1 ? active_fe_index : 0;
181 quadrature_generator.set_1D_quadrature(q1D_index);
182 quadrature_generator.generate(cell);
183
184 const Quadrature<dim> &inside_quadrature =
185 quadrature_generator.get_inside_quadrature();
186 const Quadrature<dim> &outside_quadrature =
187 quadrature_generator.get_outside_quadrature();
188 const ImmersedSurfaceQuadrature<dim> &surface_quadrature =
189 quadrature_generator.get_surface_quadrature();
190
191 // Even if a cell is formally intersected the number of created
192 // quadrature points can be 0. Avoid creating an FEValues object
193 // if that is the case.
194 if (inside_quadrature.size() > 0)
195 {
196 fe_values_inside.emplace((*mapping_collection)[mapping_index],
197 (*fe_collection)[active_fe_index],
198 inside_quadrature,
199 region_update_flags.inside);
200
201 fe_values_inside->reinit(cell);
202 }
203
204 if (outside_quadrature.size() > 0)
205 {
206 fe_values_outside.emplace((*mapping_collection)[mapping_index],
207 (*fe_collection)[active_fe_index],
208 outside_quadrature,
209 region_update_flags.outside);
210
211 fe_values_outside->reinit(cell);
212 }
213
214 if (surface_quadrature.size() > 0)
215 {
216 fe_values_surface.emplace((*mapping_collection)[mapping_index],
217 (*fe_collection)[active_fe_index],
218 surface_quadrature,
219 region_update_flags.surface);
220 fe_values_surface->reinit(cell);
221 }
222
223 break;
224 }
225 default:
226 {
228 break;
229 }
230 }
231 }
232
233
234
235 template <int dim>
236 const std::optional<::FEValues<dim>> &
238 {
239 if (current_cell_location == LocationToLevelSet::inside)
240 return fe_values_inside_full_quadrature.at(active_fe_index);
241 else
242 return fe_values_inside;
243 }
244
245
246
247 template <int dim>
248 const std::optional<::FEValues<dim>> &
250 {
251 if (current_cell_location == LocationToLevelSet::outside)
252 return fe_values_outside_full_quadrature.at(active_fe_index);
253 else
254 return fe_values_outside;
255 }
256
257
258
259 template <int dim>
260 const std::optional<FEImmersedSurfaceValues<dim>> &
262 {
263 return fe_values_surface;
264 }
265
266
267
268 template <int dim>
269 template <typename VectorType>
271 const hp::FECollection<dim> &fe_collection,
272 const Quadrature<1> &quadrature,
273 const RegionUpdateFlags region_update_flags,
274 const MeshClassifier<dim> &mesh_classifier,
275 const DoFHandler<dim> &dof_handler,
276 const VectorType &level_set,
277 const AdditionalData &additional_data)
278 : mapping_collection(&::hp::StaticMappingQ1<dim>::mapping_collection)
279 , fe_collection(&fe_collection)
280 , q_collection_1D(quadrature)
281 , region_update_flags(region_update_flags)
282 , mesh_classifier(&mesh_classifier)
283 , face_quadrature_generator(q_collection_1D,
284 dof_handler,
285 level_set,
286 additional_data)
287 {
288 // Tensor products of each quadrature in q_collection_1d. Used on the
289 // non-intersected cells.
290 hp::QCollection<dim - 1> q_collection;
291 for (const auto &quadrature : q_collection_1D)
292 q_collection.push_back(quadrature);
293
294 initialize(q_collection);
295 }
296
297
298
299 template <int dim>
300 template <typename VectorType>
302 const hp::MappingCollection<dim> &mapping_collection,
303 const hp::FECollection<dim> &fe_collection,
304 const hp::QCollection<dim - 1> &q_collection,
305 const hp::QCollection<1> &q_collection_1D,
306 const RegionUpdateFlags region_update_flags,
307 const MeshClassifier<dim> &mesh_classifier,
308 const DoFHandler<dim> &dof_handler,
309 const VectorType &level_set,
310 const AdditionalData &additional_data)
311 : mapping_collection(&mapping_collection)
312 , fe_collection(&fe_collection)
313 , q_collection_1D(q_collection_1D)
314 , region_update_flags(region_update_flags)
315 , mesh_classifier(&mesh_classifier)
316 , face_quadrature_generator(q_collection_1D,
317 dof_handler,
318 level_set,
319 additional_data)
320 {
321 initialize(q_collection);
322 }
323
324
325
326 template <int dim>
327 void
329 const hp::QCollection<dim - 1> &q_collection)
330 {
331 current_face_location = LocationToLevelSet::unassigned;
332 active_fe_index = numbers::invalid_unsigned_int;
333
334 Assert(fe_collection->size() > 0,
335 ExcMessage("Incoming hp::FECollection can not be empty."));
336 Assert(
337 mapping_collection->size() == fe_collection->size() ||
338 mapping_collection->size() == 1,
340 "Size of hp::MappingCollection must be the same as hp::FECollection or 1."));
341 Assert(
342 q_collection.size() == fe_collection->size() || q_collection.size() == 1,
344 "Size of hp::QCollection<dim> must be the same as hp::FECollection or 1."));
345 Assert(
346 q_collection_1D.size() == fe_collection->size() ||
347 q_collection_1D.size() == 1,
349 "Size of hp::QCollection<1> must be the same as hp::FECollection or 1."));
350
351 // For each element in fe_collection, create ::FEInterfaceValues
352 // objects to use on the non-intersected cells.
353 fe_values_inside_full_quadrature.resize(fe_collection->size());
354 fe_values_outside_full_quadrature.resize(fe_collection->size());
355 for (unsigned int fe_index = 0; fe_index < fe_collection->size();
356 ++fe_index)
357 {
358 const unsigned int mapping_index =
359 mapping_collection->size() > 1 ? fe_index : 0;
360 const unsigned int q_index = q_collection.size() > 1 ? fe_index : 0;
361
362 fe_values_inside_full_quadrature[fe_index].emplace(
363 (*mapping_collection)[mapping_index],
364 (*fe_collection)[fe_index],
365 q_collection[q_index],
366 region_update_flags.inside);
367 fe_values_outside_full_quadrature[fe_index].emplace(
368 (*mapping_collection)[mapping_index],
369 (*fe_collection)[fe_index],
370 q_collection[q_index],
371 region_update_flags.outside);
372 }
373 }
374
375
376
377 template <int dim>
378 template <bool level_dof_access>
379 void
382 const unsigned int face_no,
383 const std::function<void(::FEInterfaceValues<dim> &)> &call_reinit)
384 {
385 current_face_location =
386 mesh_classifier->location_to_level_set(cell, face_no);
387 active_fe_index = cell->active_fe_index();
388
389 // These objects were created with a quadrature based on the previous cell
390 // and are thus no longer valid.
391 fe_values_inside.reset();
392 fe_values_outside.reset();
393
394 switch (current_face_location)
395 {
397 {
398 call_reinit(*fe_values_inside_full_quadrature.at(active_fe_index));
399 break;
400 }
402 {
403 call_reinit(*fe_values_outside_full_quadrature.at(active_fe_index));
404 break;
405 }
407 {
408 const unsigned int mapping_index =
409 mapping_collection->size() > 1 ? active_fe_index : 0;
410 const unsigned int q1D_index =
411 q_collection_1D.size() > 1 ? active_fe_index : 0;
412
413 face_quadrature_generator.set_1D_quadrature(q1D_index);
414 face_quadrature_generator.generate(cell, face_no);
415
416 const Quadrature<dim - 1> &inside_quadrature =
417 face_quadrature_generator.get_inside_quadrature();
418 const Quadrature<dim - 1> &outside_quadrature =
419 face_quadrature_generator.get_outside_quadrature();
420
421 // Even if a cell is formally intersected the number of created
422 // quadrature points can be 0. Avoid creating an FEInterfaceValues
423 // object if that is the case.
424 if (inside_quadrature.size() > 0)
425 {
426 fe_values_inside.emplace((*mapping_collection)[mapping_index],
427 (*fe_collection)[active_fe_index],
428 inside_quadrature,
429 region_update_flags.inside);
430
431 call_reinit(*fe_values_inside);
432 }
433
434 if (outside_quadrature.size() > 0)
435 {
436 fe_values_outside.emplace((*mapping_collection)[mapping_index],
437 (*fe_collection)[active_fe_index],
438 outside_quadrature,
439 region_update_flags.outside);
440
441 call_reinit(*fe_values_outside);
442 }
443 break;
444 }
445 default:
446 {
448 break;
449 }
450 }
451 }
452
453
454
455 template <int dim>
456 const std::optional<::FEInterfaceValues<dim>> &
458 {
459 if (current_face_location == LocationToLevelSet::inside)
460 return fe_values_inside_full_quadrature.at(active_fe_index);
461 else
462 return fe_values_inside;
463 }
464
465
466
467 template <int dim>
468 const std::optional<::FEInterfaceValues<dim>> &
470 {
471 if (current_face_location == LocationToLevelSet::outside)
472 return fe_values_outside_full_quadrature.at(active_fe_index);
473 else
474 return fe_values_outside;
475 }
476
477
478#include "fe_values.inst"
479
480} // namespace NonMatching
const std::optional<::FEInterfaceValues< dim > > & get_outside_fe_values() const
Definition fe_values.cc:469
void initialize(const hp::QCollection< dim - 1 > &q_collection)
Definition fe_values.cc:328
const hp::QCollection< 1 > q_collection_1D
Definition fe_values.h:602
typename FaceQuadratureGenerator< dim >::AdditionalData AdditionalData
Definition fe_values.h:439
const std::optional<::FEInterfaceValues< dim > > & get_inside_fe_values() const
Definition fe_values.cc:457
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:270
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:380
const std::optional<::FEValues< dim > > & get_outside_fe_values() const
Definition fe_values.cc:249
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:261
const hp::QCollection< 1 > q_collection_1D
Definition fe_values.h:283
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access > > &cell)
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:237
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
#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