Reference documentation for deal.II version GIT relicensing-233-g802318d791 2024-03-28 20:20: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
loop.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2009 - 2023 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
15
16#ifndef dealii_mesh_worker_loop_h
17#define dealii_mesh_worker_loop_h
18
19#include <deal.II/base/config.h>
20
23
25#include <deal.II/grid/tria.h>
26
30
31#include <functional>
32
34
35// Forward declaration
36#ifndef DOXYGEN
37template <typename>
39#endif
40
41namespace internal
42{
46 template <class DI>
47 inline bool
49 {
50 return false;
51 }
52
53 template <typename AccessorType>
54 inline bool
56 {
57 return true;
58 }
59
60 template <typename AccessorType>
61 inline bool
63 const ::FilteredIterator<TriaActiveIterator<AccessorType>> &)
64 {
65 return true;
66 }
67
68 template <int dim, class DOFINFO, class A>
69 void
71 {
72 dinfo.assemble(*assembler);
73 }
74} // namespace internal
75
76
77
78namespace MeshWorker
79{
84 {
85 public:
90 : own_cells(true)
91 , ghost_cells(false)
94 , cells_first(true)
95 {}
96
101
107
128
142
153
159 };
160
161
162
188 template <class INFOBOX,
189 class DOFINFO,
190 int dim,
191 int spacedim,
192 typename IteratorType>
193 void
195 IteratorType cell,
196 DoFInfoBox<dim, DOFINFO> &dof_info,
197 INFOBOX &info,
198 const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
199 &cell_worker,
200 const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
201 &boundary_worker,
202 const std::function<void(DOFINFO &,
203 DOFINFO &,
204 typename INFOBOX::CellInfo &,
205 typename INFOBOX::CellInfo &)> &face_worker,
206 const LoopControl &loop_control)
207 {
208 const bool ignore_subdomain =
209 (cell->get_triangulation().locally_owned_subdomain() ==
211
212 types::subdomain_id csid = (cell->is_level_cell()) ?
213 cell->level_subdomain_id() :
214 cell->subdomain_id();
215
216 const bool own_cell =
217 ignore_subdomain ||
218 (csid == cell->get_triangulation().locally_owned_subdomain());
219
220 dof_info.reset();
221
222 if ((!ignore_subdomain) && (csid == numbers::artificial_subdomain_id))
223 return;
224
225 dof_info.cell.reinit(cell);
226 dof_info.cell_valid = true;
227
228 const bool integrate_cell = (cell_worker != nullptr);
229 const bool integrate_boundary = (boundary_worker != nullptr);
230 const bool integrate_interior_face = (face_worker != nullptr);
231
232 if (integrate_cell)
233 info.cell.reinit(dof_info.cell);
234 // Execute this, if cells
235 // have to be dealt with
236 // before faces
237 if (integrate_cell && loop_control.cells_first &&
238 ((loop_control.own_cells && own_cell) ||
239 (loop_control.ghost_cells && !own_cell)))
240 cell_worker(dof_info.cell, info.cell);
241
242 // Call the callback function in
243 // the info box to do
244 // computations between cell and
245 // face action.
246 info.post_cell(dof_info);
247
248 if (integrate_interior_face || integrate_boundary)
249 for (const unsigned int face_no : cell->face_indices())
250 {
251 typename IteratorType::AccessorType::Container::face_iterator face =
252 cell->face(face_no);
253 if (cell->at_boundary(face_no) &&
254 !cell->has_periodic_neighbor(face_no))
255 {
256 // only integrate boundary faces of own cells
257 if (integrate_boundary && own_cell)
258 {
259 dof_info.interior_face_available[face_no] = true;
260 dof_info.interior[face_no].reinit(cell, face, face_no);
261 info.boundary.reinit(dof_info.interior[face_no]);
262 boundary_worker(dof_info.interior[face_no], info.boundary);
263 }
264 }
265 else if (integrate_interior_face)
266 {
267 // Interior face
269 cell->neighbor_or_periodic_neighbor(face_no);
270
272 if (neighbor->is_level_cell())
273 neighbid = neighbor->level_subdomain_id();
274 // subdomain id is only valid for active cells
275 else if (neighbor->is_active())
276 neighbid = neighbor->subdomain_id();
277
278 const bool own_neighbor =
279 ignore_subdomain ||
280 (neighbid ==
281 cell->get_triangulation().locally_owned_subdomain());
282
283 // skip all faces between two ghost cells
284 if (!own_cell && !own_neighbor)
285 continue;
286
287 // skip if the user doesn't want faces between own cells
288 if (own_cell && own_neighbor &&
289 loop_control.own_faces == LoopControl::never)
290 continue;
291
292 // skip face to ghost
293 if (own_cell != own_neighbor &&
294 loop_control.faces_to_ghost == LoopControl::never)
295 continue;
296
297 // Deal with refinement edges from the refined side. Assuming
298 // one-irregular meshes, this situation should only occur if both
299 // cells are active.
300 const bool periodic_neighbor =
301 cell->has_periodic_neighbor(face_no);
302
303 if ((!periodic_neighbor && cell->neighbor_is_coarser(face_no)) ||
304 (periodic_neighbor &&
305 cell->periodic_neighbor_is_coarser(face_no)))
306 {
307 Assert(cell->is_active(), ExcInternalError());
308 Assert(neighbor->is_active(), ExcInternalError());
309
310 // skip if only one processor needs to assemble the face
311 // to a ghost cell and the fine cell is not ours.
312 if (!own_cell &&
313 loop_control.faces_to_ghost == LoopControl::one)
314 continue;
315
316 const std::pair<unsigned int, unsigned int> neighbor_face_no =
317 periodic_neighbor ?
318 cell->periodic_neighbor_of_coarser_periodic_neighbor(
319 face_no) :
320 cell->neighbor_of_coarser_neighbor(face_no);
321 const typename IteratorType::AccessorType::Container::
322 face_iterator nface =
323 neighbor->face(neighbor_face_no.first);
324
325 dof_info.interior_face_available[face_no] = true;
326 dof_info.exterior_face_available[face_no] = true;
327 dof_info.interior[face_no].reinit(cell, face, face_no);
328 info.face.reinit(dof_info.interior[face_no]);
329 dof_info.exterior[face_no].reinit(neighbor,
330 nface,
331 neighbor_face_no.first,
332 neighbor_face_no.second);
333 info.subface.reinit(dof_info.exterior[face_no]);
334
335 face_worker(dof_info.interior[face_no],
336 dof_info.exterior[face_no],
337 info.face,
338 info.subface);
339 }
340 else
341 {
342 // If iterator is active and neighbor is refined, skip
343 // internal face.
345 neighbor->has_children())
346 {
347 Assert(
348 loop_control.own_faces != LoopControl::both,
350 "Assembling from both sides for own_faces is not "
351 "supported with hanging nodes!"));
352 continue;
353 }
354
355 // Now neighbor is on same level, double-check this:
356 Assert(cell->level() == neighbor->level(),
358
359 // If we own both cells only do faces from one side (unless
360 // LoopControl says otherwise). Here, we rely on cell
361 // comparison that will look at cell->index().
362 if (own_cell && own_neighbor &&
363 loop_control.own_faces == LoopControl::one &&
364 (neighbor < cell))
365 continue;
366
367 // independent of loop_control.faces_to_ghost,
368 // we only look at faces to ghost on the same level once
369 // (only where own_cell=true and own_neighbor=false)
370 if (!own_cell)
371 continue;
372
373 // now only one processor assembles faces_to_ghost. We let the
374 // processor with the smaller (level-)subdomain id assemble
375 // the face.
376 if (own_cell && !own_neighbor &&
377 loop_control.faces_to_ghost == LoopControl::one &&
378 (neighbid < csid))
379 continue;
380
381 const unsigned int neighbor_face_no =
382 periodic_neighbor ?
383 cell->periodic_neighbor_face_no(face_no) :
384 cell->neighbor_face_no(face_no);
385 Assert(periodic_neighbor ||
386 neighbor->face(neighbor_face_no) == face,
388 // Regular interior face
389 dof_info.interior_face_available[face_no] = true;
390 dof_info.exterior_face_available[face_no] = true;
391 dof_info.interior[face_no].reinit(cell, face, face_no);
392 info.face.reinit(dof_info.interior[face_no]);
393 dof_info.exterior[face_no].reinit(neighbor,
394 neighbor->face(
395 neighbor_face_no),
396 neighbor_face_no);
397 info.neighbor.reinit(dof_info.exterior[face_no]);
398
399 face_worker(dof_info.interior[face_no],
400 dof_info.exterior[face_no],
401 info.face,
402 info.neighbor);
403 }
404 }
405 } // faces
406 // Call the callback function in
407 // the info box to do
408 // computations between face and
409 // cell action.
410 info.post_faces(dof_info);
411
412 // Execute this, if faces
413 // have to be handled first
414 if (integrate_cell && !loop_control.cells_first &&
415 ((loop_control.own_cells && own_cell) ||
416 (loop_control.ghost_cells && !own_cell)))
417 cell_worker(dof_info.cell, info.cell);
418 }
419
420
435 template <int dim,
436 int spacedim,
437 class DOFINFO,
438 class INFOBOX,
439 typename AssemblerType,
440 typename IteratorType>
441 void
442 loop(IteratorType begin,
444 DOFINFO &dinfo,
445 INFOBOX &info,
446 const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
447 &cell_worker,
448 const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
449 &boundary_worker,
450 const std::function<void(DOFINFO &,
451 DOFINFO &,
452 typename INFOBOX::CellInfo &,
453 typename INFOBOX::CellInfo &)> &face_worker,
454 AssemblerType &assembler,
455 const LoopControl &lctrl = LoopControl())
456 {
457 DoFInfoBox<dim, DOFINFO> dof_info(dinfo);
458
459 assembler.initialize_info(dof_info.cell, false);
460 for (const unsigned int i : GeometryInfo<dim>::face_indices())
461 {
462 assembler.initialize_info(dof_info.interior[i], true);
463 assembler.initialize_info(dof_info.exterior[i], true);
464 }
465
466 // Loop over all cells
468 begin,
469 end,
470 [&cell_worker, &boundary_worker, &face_worker, &lctrl](
471 IteratorType cell, INFOBOX &info, DoFInfoBox<dim, DOFINFO> &dof_info) {
472 cell_action<INFOBOX, DOFINFO, dim, spacedim, IteratorType>(
473 cell,
474 dof_info,
475 info,
476 cell_worker,
477 boundary_worker,
478 face_worker,
479 lctrl);
480 },
481 [&assembler](const MeshWorker::DoFInfoBox<dim, DOFINFO> &dinfo) {
482 ::internal::assemble<dim, DOFINFO, AssemblerType>(dinfo,
483 &assembler);
484 },
485 info,
486 dof_info);
487 }
488
489
496 template <int dim,
497 int spacedim,
498 typename IteratorType,
499 typename AssemblerType>
500 void
501 integration_loop(IteratorType begin,
503 DoFInfo<dim, spacedim> &dof_info,
505 const LocalIntegrator<dim, spacedim> &integrator,
506 AssemblerType &assembler,
507 const LoopControl &lctrl = LoopControl())
508 {
509 std::function<void(DoFInfo<dim, spacedim> &,
511 cell_worker;
512 std::function<void(DoFInfo<dim, spacedim> &,
514 boundary_worker;
515 std::function<void(DoFInfo<dim, spacedim> &,
519 face_worker;
520 if (integrator.use_cell)
521 cell_worker =
522 [&integrator](DoFInfo<dim, spacedim> &dof_info,
523 IntegrationInfo<dim, spacedim> &integration_info) {
524 integrator.cell(dof_info, integration_info);
525 };
526 if (integrator.use_boundary)
527 boundary_worker =
528 [&integrator](DoFInfo<dim, spacedim> &dof_info,
529 IntegrationInfo<dim, spacedim> &integration_info) {
530 integrator.boundary(dof_info, integration_info);
531 };
532 if (integrator.use_face)
533 face_worker =
534 [&integrator](DoFInfo<dim, spacedim> &dof_info_1,
535 DoFInfo<dim, spacedim> &dof_info_2,
536 IntegrationInfo<dim, spacedim> &integration_info_1,
537 IntegrationInfo<dim, spacedim> &integration_info_2) {
538 integrator.face(dof_info_1,
539 dof_info_2,
540 integration_info_1,
541 integration_info_2);
542 };
543
544 loop<dim, spacedim>(begin,
545 end,
546 dof_info,
547 box,
548 cell_worker,
549 boundary_worker,
550 face_worker,
551 assembler,
552 lctrl);
553 }
554
555} // namespace MeshWorker
556
558
559#endif
DOFINFO interior[GeometryInfo< dim >::faces_per_cell]
Definition dof_info.h:263
void assemble(ASSEMBLER &ass) const
Definition dof_info.h:492
bool exterior_face_available[GeometryInfo< dim >::faces_per_cell]
Definition dof_info.h:279
bool interior_face_available[GeometryInfo< dim >::faces_per_cell]
Definition dof_info.h:273
DOFINFO exterior[GeometryInfo< dim >::faces_per_cell]
Definition dof_info.h:267
virtual void cell(DoFInfo< dim, spacedim, number > &dinfo, IntegrationInfo< dim, spacedim > &info) const
virtual void face(DoFInfo< dim, spacedim, number > &dinfo1, DoFInfo< dim, spacedim, number > &dinfo2, IntegrationInfo< dim, spacedim > &info1, IntegrationInfo< dim, spacedim > &info2) const
virtual void boundary(DoFInfo< dim, spacedim, number > &dinfo, IntegrationInfo< dim, spacedim > &info) const
FaceOption own_faces
Definition loop.h:152
FaceOption faces_to_ghost
Definition loop.h:141
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
void cell_action(IteratorType cell, DoFInfoBox< dim, DOFINFO > &dof_info, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, const LoopControl &loop_control)
Definition loop.h:194
void integration_loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DoFInfo< dim, spacedim > &dof_info, IntegrationInfoBox< dim, spacedim > &box, const LocalIntegrator< dim, spacedim > &integrator, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:501
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:442
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
bool is_active_iterator(const DI &)
Definition loop.h:48
const types::subdomain_id artificial_subdomain_id
Definition types.h:362
const types::subdomain_id invalid_subdomain_id
Definition types.h:341
typename type_identity< T >::type type_identity_t
Definition type_traits.h:95
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()