Reference documentation for deal.II version 9.5.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// 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
18
29#include <deal.II/lac/vector.h>
30
32
34
36
37namespace 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 (const auto &quadrature : q_collection_1D)
70 q_collection.push_back(Quadrature<dim>(quadrature));
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 (const auto &quadrature : q_collection_1D)
296 q_collection.push_back(quadrature);
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,
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,
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,
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
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 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
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
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access > > &cell)
Definition fe_values.cc:154
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:265
void push_back(const Quadrature< dim_in > &new_quadrature)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
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