Reference documentation for deal.II version GIT relicensing-891-g13379fae33 2024-06-22 22: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
tria.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 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
15
18#include <deal.II/base/mpi.templates.h>
23
25
30#include <deal.II/grid/tria.h>
35
36#include <boost/archive/text_iarchive.hpp>
37#include <boost/archive/text_oarchive.hpp>
38
39#include <algorithm>
40#include <array>
41#include <cmath>
42#include <cstdint>
43#include <fstream>
44#include <functional>
45#include <limits>
46#include <list>
47#include <map>
48#include <memory>
49#include <numeric>
50
51
53
54
55namespace internal
56{
57 namespace TriangulationImplementation
58 {
60 : n_levels(0)
61 , n_lines(0)
62 , n_active_lines(0)
63 // all other fields are
64 // default constructed
65 {}
66
67
68
69 std::size_t
71 {
72 std::size_t mem =
77 MemoryConsumption::memory_consumption(n_active_lines_level);
78
79 if (active_cell_index_partitioner)
80 mem += active_cell_index_partitioner->memory_consumption();
81
82 for (const auto &partitioner : level_cell_index_partitioners)
83 if (partitioner)
84 mem += partitioner->memory_consumption();
85
86 return mem;
87 }
88
89
91 : n_quads(0)
92 , n_active_quads(0)
93 // all other fields are
94 // default constructed
95 {}
96
97
98
99 std::size_t
108
109
110
112 : n_hexes(0)
113 , n_active_hexes(0)
114 // all other fields are
115 // default constructed
116 {}
117
118
119
120 std::size_t
129 } // namespace TriangulationImplementation
130
131
132 template <int dim, int spacedim>
135 : variable_size_data_stored(false)
136 {}
137
138
139 template <int dim, int spacedim>
141 void CellAttachedDataSerializer<dim, spacedim>::pack_data(
142 const std::vector<cell_relation_t> &cell_relations,
143 const std::vector<
144 typename internal::CellAttachedData<dim, spacedim>::pack_callback_t>
145 &pack_callbacks_fixed,
146 const std::vector<
147 typename internal::CellAttachedData<dim, spacedim>::pack_callback_t>
148 &pack_callbacks_variable,
149 const MPI_Comm &mpi_communicator)
150 {
151 Assert(src_data_fixed.empty(),
152 ExcMessage("Previously packed data has not been released yet!"));
153 Assert(src_sizes_variable.empty(), ExcInternalError());
154
155 const unsigned int n_callbacks_fixed = pack_callbacks_fixed.size();
156 const unsigned int n_callbacks_variable = pack_callbacks_variable.size();
157
158 // Store information that we packed variable size data in
159 // a member variable for later.
160 variable_size_data_stored = (n_callbacks_variable > 0);
161
162 // If variable transfer is scheduled, we will store the data size that
163 // each variable size callback function writes in this auxiliary
164 // container. The information will be stored by each cell in this vector
165 // temporarily.
166 std::vector<unsigned int> cell_sizes_variable_cumulative(
167 n_callbacks_variable);
168
169 // Prepare the buffer structure, in which each callback function will
170 // store its data for each active cell.
171 // The outmost shell in this container construct corresponds to the
172 // data packed per cell. The next layer resembles the data that
173 // each callback function packs on the corresponding cell. These
174 // buffers are chains of chars stored in an std::vector<char>.
175 // A visualisation of the data structure:
176 /* clang-format off */
177 // | cell_1 | | cell_2 | ...
178 // || callback_1 || callback_2 |...| || callback_1 || callback_2 |...| ...
179 // |||char|char|...|||char|char|...|...| |||char|char|...|||char|char|...|...| ...
180 /* clang-format on */
181 std::vector<std::vector<std::vector<char>>> packed_fixed_size_data(
182 cell_relations.size());
183 std::vector<std::vector<std::vector<char>>> packed_variable_size_data(
184 variable_size_data_stored ? cell_relations.size() : 0);
185
186 //
187 // --------- Pack data for fixed and variable size transfer ---------
188 //
189 // Iterate over all cells, call all callback functions on each cell,
190 // and store their data in the corresponding buffer scope.
191 {
192 auto cell_rel_it = cell_relations.cbegin();
193 auto data_cell_fixed_it = packed_fixed_size_data.begin();
194 auto data_cell_variable_it = packed_variable_size_data.begin();
195 for (; cell_rel_it != cell_relations.cend(); ++cell_rel_it)
196 {
197 const auto &dealii_cell = cell_rel_it->first;
198 const auto &cell_status = cell_rel_it->second;
199
200 // Assertions about the tree structure.
201 switch (cell_status)
202 {
205 // double check the condition that we will only ever attach
206 // data to active cells when we get here
207 Assert(dealii_cell->is_active(), ExcInternalError());
208 break;
209
211 // double check the condition that we will only ever attach
212 // data to cells with children when we get here. however, we
213 // can only tolerate one level of coarsening at a time, so
214 // check that the children are all active
215 Assert(dealii_cell->is_active() == false, ExcInternalError());
216 for (unsigned int c = 0;
217 c < GeometryInfo<dim>::max_children_per_cell;
218 ++c)
219 Assert(dealii_cell->child(c)->is_active(),
221 break;
222
224 // do nothing on invalid cells
225 break;
226
227 default:
229 break;
230 }
231
232 // Reserve memory corresponding to the number of callback
233 // functions that will be called.
234 // If variable size transfer is scheduled, we need to leave
235 // room for an array that holds information about how many
236 // bytes each of the variable size callback functions will
237 // write.
238 // On cells flagged with CellStatus::cell_invalid, only its CellStatus
239 // will be stored.
240 const unsigned int n_fixed_size_data_sets_on_cell =
241 1 + ((cell_status == CellStatus::cell_invalid) ?
242 0 :
243 ((variable_size_data_stored ? 1 : 0) + n_callbacks_fixed));
244 data_cell_fixed_it->resize(n_fixed_size_data_sets_on_cell);
245
246 // We continue with packing all data on this specific cell.
247 auto data_fixed_it = data_cell_fixed_it->begin();
248
249 // First, we pack the CellStatus information.
250 // to get consistent data sizes on each cell for the fixed size
251 // transfer, we won't allow compression
252 *data_fixed_it =
253 Utilities::pack(cell_status, /*allow_compression=*/false);
254 ++data_fixed_it;
255
256 // Proceed with all registered callback functions.
257 // Skip cells with the CellStatus::cell_invalid flag.
258 if (cell_status != CellStatus::cell_invalid)
259 {
260 // Pack fixed size data.
261 for (auto callback_it = pack_callbacks_fixed.cbegin();
262 callback_it != pack_callbacks_fixed.cend();
263 ++callback_it, ++data_fixed_it)
264 {
265 *data_fixed_it = (*callback_it)(dealii_cell, cell_status);
266 }
267
268 // Pack variable size data.
269 // If we store variable size data, we need to transfer
270 // the sizes of each corresponding callback function
271 // via fixed size transfer as well.
272 if (variable_size_data_stored)
273 {
274 const unsigned int n_variable_size_data_sets_on_cell =
275 ((cell_status == CellStatus::cell_invalid) ?
276 0 :
277 n_callbacks_variable);
278 data_cell_variable_it->resize(
279 n_variable_size_data_sets_on_cell);
280
281 auto callback_it = pack_callbacks_variable.cbegin();
282 auto data_variable_it = data_cell_variable_it->begin();
283 auto sizes_variable_it =
284 cell_sizes_variable_cumulative.begin();
285 for (; callback_it != pack_callbacks_variable.cend();
286 ++callback_it, ++data_variable_it, ++sizes_variable_it)
287 {
288 *data_variable_it =
289 (*callback_it)(dealii_cell, cell_status);
290
291 // Store data sizes for each callback function first.
292 // Make it cumulative below.
293 *sizes_variable_it = data_variable_it->size();
294 }
295
296 // Turn size vector into its cumulative representation.
297 std::partial_sum(cell_sizes_variable_cumulative.begin(),
298 cell_sizes_variable_cumulative.end(),
299 cell_sizes_variable_cumulative.begin());
300
301 // Serialize cumulative variable size vector
302 // value-by-value. This way we can circumvent the overhead
303 // of storing the container object as a whole, since we
304 // know its size by the number of registered callback
305 // functions.
306 data_fixed_it->resize(n_callbacks_variable *
307 sizeof(unsigned int));
308 for (unsigned int i = 0; i < n_callbacks_variable; ++i)
309 std::memcpy(&(data_fixed_it->at(i * sizeof(unsigned int))),
310 &(cell_sizes_variable_cumulative.at(i)),
311 sizeof(unsigned int));
312
313 ++data_fixed_it;
314 }
315
316 // Double check that we packed everything we wanted
317 // in the fixed size buffers.
318 Assert(data_fixed_it == data_cell_fixed_it->end(),
320 }
321
322 ++data_cell_fixed_it;
323
324 // Increment the variable size data iterator
325 // only if we actually pack this kind of data
326 // to avoid getting out of bounds.
327 if (variable_size_data_stored)
328 ++data_cell_variable_it;
329 } // loop over cell_relations
330 }
331
332 //
333 // ----------- Gather data sizes for fixed size transfer ------------
334 //
335 // Generate a vector which stores the sizes of each callback function,
336 // including the packed CellStatus transfer.
337 // Find the very first cell that we wrote to with all callback
338 // functions (i.e. a cell that was not flagged with
339 // CellStatus::cell_invalid) and store the sizes of each buffer.
340 //
341 // To deal with the case that at least one of the processors does not
342 // own any cell at all, we will exchange the information about the data
343 // sizes among them later. The code in between is still well-defined,
344 // since the following loops will be skipped.
345 std::vector<unsigned int> local_sizes_fixed(
346 1 + n_callbacks_fixed + (variable_size_data_stored ? 1 : 0));
347 for (const auto &data_cell : packed_fixed_size_data)
348 {
349 if (data_cell.size() == local_sizes_fixed.size())
350 {
351 auto sizes_fixed_it = local_sizes_fixed.begin();
352 auto data_fixed_it = data_cell.cbegin();
353 for (; data_fixed_it != data_cell.cend();
354 ++data_fixed_it, ++sizes_fixed_it)
355 {
356 *sizes_fixed_it = data_fixed_it->size();
357 }
358
359 break;
360 }
361 }
362
363 // Check if all cells have valid sizes.
364 for (auto data_cell_fixed_it = packed_fixed_size_data.cbegin();
365 data_cell_fixed_it != packed_fixed_size_data.cend();
366 ++data_cell_fixed_it)
367 {
368 Assert((data_cell_fixed_it->size() == 1) ||
369 (data_cell_fixed_it->size() == local_sizes_fixed.size()),
371 }
372
373 // Share information about the packed data sizes
374 // of all callback functions across all processors, in case one
375 // of them does not own any cells at all.
376 std::vector<unsigned int> global_sizes_fixed(local_sizes_fixed.size());
377 Utilities::MPI::max(local_sizes_fixed,
378 mpi_communicator,
379 global_sizes_fixed);
380
381 // Construct cumulative sizes, since this is the only information
382 // we need from now on.
383 sizes_fixed_cumulative.resize(global_sizes_fixed.size());
384 std::partial_sum(global_sizes_fixed.begin(),
385 global_sizes_fixed.end(),
386 sizes_fixed_cumulative.begin());
387
388 //
389 // ---------- Gather data sizes for variable size transfer ----------
390 //
391 if (variable_size_data_stored)
392 {
393 src_sizes_variable.reserve(packed_variable_size_data.size());
394 for (const auto &data_cell : packed_variable_size_data)
395 {
396 int variable_data_size_on_cell = 0;
397
398 for (const auto &data : data_cell)
399 variable_data_size_on_cell += data.size();
401 src_sizes_variable.push_back(variable_data_size_on_cell);
402 }
403 }
404
405 //
406 // ------------------------ Build buffers ---------------------------
407 //
408 const unsigned int expected_size_fixed =
409 cell_relations.size() * sizes_fixed_cumulative.back();
410 const unsigned int expected_size_variable =
411 std::accumulate(src_sizes_variable.begin(),
412 src_sizes_variable.end(),
413 std::vector<int>::size_type(0));
414
415 // Move every piece of packed fixed size data into the consecutive
416 // buffer.
417 src_data_fixed.reserve(expected_size_fixed);
418 for (const auto &data_cell_fixed : packed_fixed_size_data)
419 {
420 // Move every fraction of packed data into the buffer
421 // reserved for this particular cell.
422 for (const auto &data_fixed : data_cell_fixed)
423 std::move(data_fixed.begin(),
424 data_fixed.end(),
425 std::back_inserter(src_data_fixed));
426
427 // If we only packed the CellStatus information
428 // (i.e. encountered a cell flagged CellStatus::cell_invalid),
429 // fill the remaining space with invalid entries.
430 // We can skip this if there is nothing else to pack.
431 if ((data_cell_fixed.size() == 1) &&
432 (sizes_fixed_cumulative.size() > 1))
434 const std::size_t bytes_skipped =
435 sizes_fixed_cumulative.back() - sizes_fixed_cumulative.front();
436
437 src_data_fixed.insert(src_data_fixed.end(),
438 bytes_skipped,
439 static_cast<char>(-1)); // invalid_char
440 }
441 }
442
443 // Move every piece of packed variable size data into the consecutive
444 // buffer.
445 if (variable_size_data_stored)
446 {
447 src_data_variable.reserve(expected_size_variable);
448 for (const auto &data_cell : packed_variable_size_data)
449 {
450 // Move every fraction of packed data into the buffer
451 // reserved for this particular cell.
452 for (const auto &data : data_cell)
453 std::move(data.begin(),
454 data.end(),
455 std::back_inserter(src_data_variable));
456 }
458
459 // Double check that we packed everything correctly.
460 Assert(src_data_fixed.size() == expected_size_fixed, ExcInternalError());
461 Assert(src_data_variable.size() == expected_size_variable,
463 }
464
465
466
467 template <int dim, int spacedim>
469 void CellAttachedDataSerializer<dim, spacedim>::unpack_cell_status(
470 std::vector<
471 typename CellAttachedDataSerializer<dim, spacedim>::cell_relation_t>
472 &cell_relations) const
473 {
474 Assert(sizes_fixed_cumulative.size() > 0,
475 ExcMessage("No data has been packed!"));
476 if (cell_relations.size() > 0)
477 {
478 Assert(dest_data_fixed.size() > 0,
479 ExcMessage("No data has been received!"));
480 }
482 // Size of CellStatus object that will be unpacked on each cell.
483 const unsigned int size = sizes_fixed_cumulative.front();
484
485 // Iterate over all cells and overwrite the CellStatus
486 // information from the transferred data.
487 // Proceed buffer iterator position to next cell after
488 // each iteration.
489 auto cell_rel_it = cell_relations.begin();
490 auto dest_fixed_it = dest_data_fixed.cbegin();
491 for (; cell_rel_it != cell_relations.end();
492 ++cell_rel_it, dest_fixed_it += sizes_fixed_cumulative.back())
493 {
494 cell_rel_it->second = // cell_status
495 Utilities::unpack<CellStatus>(dest_fixed_it,
496 dest_fixed_it + size,
497 /*allow_compression=*/false);
498 }
499 }
500
501
502
503 template <int dim, int spacedim>
505 void CellAttachedDataSerializer<dim, spacedim>::unpack_data(
506 const std::vector<
507 typename CellAttachedDataSerializer<dim, spacedim>::cell_relation_t>
508 &cell_relations,
509 const unsigned int handle,
510 const std::function<
511 void(const cell_iterator &,
512 const CellStatus &,
513 const boost::iterator_range<std::vector<char>::const_iterator> &)>
514 &unpack_callback) const
515 {
516 // We decode the handle returned by register_data_attach() back into
517 // a format we can use. All even handles belong to those callback
518 // functions which write/read variable size data, all odd handles
519 // interact with fixed size buffers.
520 const bool callback_variable_transfer = (handle % 2 == 0);
521 const unsigned int callback_index = handle / 2;
522
523 // Cells will always receive fixed size data (i.e., CellStatus
524 // information), but not necessarily variable size data (e.g., with a
525 // ParticleHandler a cell might not contain any particle at all).
526 // Thus it is sufficient to check if fixed size data has been received.
527 Assert(sizes_fixed_cumulative.size() > 0,
528 ExcMessage("No data has been packed!"));
529 if (cell_relations.size() > 0)
530 {
531 Assert(dest_data_fixed.size() > 0,
532 ExcMessage("No data has been received!"));
533 }
534
535 std::vector<char>::const_iterator dest_data_it;
536 std::vector<char>::const_iterator dest_sizes_cell_it;
537
538 // Depending on whether our callback function unpacks fixed or
539 // variable size data, we have to pursue different approaches
540 // to localize the correct fraction of the buffer from which
541 // we are allowed to read.
542 unsigned int offset = numbers::invalid_unsigned_int;
543 unsigned int size = numbers::invalid_unsigned_int;
544 unsigned int data_increment = numbers::invalid_unsigned_int;
545
546 if (callback_variable_transfer)
547 {
548 // For the variable size data, we need to extract the
549 // data size from the fixed size buffer on each cell.
550 //
551 // We packed this information last, so the last packed
552 // object in the fixed size buffer corresponds to the
553 // variable data sizes.
554 //
555 // The last entry of sizes_fixed_cumulative corresponds
556 // to the size of all fixed size data packed on the cell.
557 // To get the offset for the last packed object, we need
558 // to get the next-to-last entry.
559 const unsigned int offset_variable_data_sizes =
560 sizes_fixed_cumulative[sizes_fixed_cumulative.size() - 2];
561
562 // This iterator points to the data size that the
563 // callback_function packed for each specific cell.
564 // Adjust buffer iterator to the offset of the callback
565 // function so that we only have to advance its position
566 // to the next cell after each iteration.
567 dest_sizes_cell_it = dest_data_fixed.cbegin() +
568 offset_variable_data_sizes +
569 callback_index * sizeof(unsigned int);
570
571 // Let the data iterator point to the correct buffer.
572 dest_data_it = dest_data_variable.cbegin();
573 }
574 else
575 {
576 // For the fixed size data, we can get the information about
577 // the buffer location on each cell directly from the
578 // sizes_fixed_cumulative vector.
579 offset = sizes_fixed_cumulative[callback_index];
580 size = sizes_fixed_cumulative[callback_index + 1] - offset;
581 data_increment = sizes_fixed_cumulative.back();
582
583 // Let the data iterator point to the correct buffer.
584 // Adjust buffer iterator to the offset of the callback
585 // function so that we only have to advance its position
586 // to the next cell after each iteration.
587 if (cell_relations.begin() != cell_relations.end())
588 dest_data_it = dest_data_fixed.cbegin() + offset;
589 }
590
591 // Iterate over all cells and unpack the transferred data.
592 auto cell_rel_it = cell_relations.begin();
593 auto dest_sizes_it = dest_sizes_variable.cbegin();
594 for (; cell_rel_it != cell_relations.end(); ++cell_rel_it)
595 {
596 const auto &dealii_cell = cell_rel_it->first;
597 const auto &cell_status = cell_rel_it->second;
598
599 if (callback_variable_transfer)
600 {
601 // Update the increment according to the whole data size
602 // of the current cell.
603 data_increment = *dest_sizes_it;
604
605 if (cell_status != CellStatus::cell_invalid)
606 {
607 // Extract the corresponding values for offset and size from
608 // the cumulative sizes array stored in the fixed size
609 // buffer.
610 if (callback_index == 0)
611 offset = 0;
612 else
613 std::memcpy(&offset,
614 &(*(dest_sizes_cell_it - sizeof(unsigned int))),
615 sizeof(unsigned int));
616
617 std::memcpy(&size,
618 &(*dest_sizes_cell_it),
619 sizeof(unsigned int));
620
621 size -= offset;
622
623 // Move the data iterator to the corresponding position
624 // of the callback function and adjust the increment
625 // accordingly.
626 dest_data_it += offset;
627 data_increment -= offset;
628 }
629
630 // Advance data size iterators to the next cell, avoid iterating
631 // past the end of dest_sizes_cell_it
632 if (cell_rel_it != cell_relations.end() - 1)
633 dest_sizes_cell_it += sizes_fixed_cumulative.back();
634 ++dest_sizes_it;
635 }
636
637 switch (cell_status)
638 {
641 unpack_callback(dealii_cell,
642 cell_status,
643 boost::make_iterator_range(dest_data_it,
644 dest_data_it + size));
645 break;
646
648 unpack_callback(dealii_cell->parent(),
649 cell_status,
650 boost::make_iterator_range(dest_data_it,
651 dest_data_it + size));
652 break;
653
655 // Skip this cell.
656 break;
657
658 default:
660 break;
661 }
662
663 if (cell_rel_it != cell_relations.end() - 1)
664 dest_data_it += data_increment;
665 }
666 }
667
668
669
670 template <int dim, int spacedim>
672 void CellAttachedDataSerializer<dim, spacedim>::save(
673 const unsigned int global_first_cell,
674 const unsigned int global_num_cells,
675 const std::string &filename,
676 const MPI_Comm &mpi_communicator) const
677 {
678 Assert(sizes_fixed_cumulative.size() > 0,
679 ExcMessage("No data has been packed!"));
680
681#ifdef DEAL_II_WITH_MPI
682 // Large fractions of this function have been copied from
683 // DataOutInterface::write_vtu_in_parallel.
684 // TODO: Write general MPIIO interface.
685
686 const int myrank = Utilities::MPI::this_mpi_process(mpi_communicator);
687 const int mpisize = Utilities::MPI::n_mpi_processes(mpi_communicator);
688
689 if (mpisize > 1)
690 {
691 const unsigned int bytes_per_cell = sizes_fixed_cumulative.back();
692
693 //
694 // ---------- Fixed size data ----------
695 //
696 {
697 const std::string fname_fixed = std::string(filename) + "_fixed.data";
698
699 MPI_Info info;
700 int ierr = MPI_Info_create(&info);
701 AssertThrowMPI(ierr);
702
703 MPI_File fh;
704 ierr = MPI_File_open(mpi_communicator,
705 fname_fixed.c_str(),
706 MPI_MODE_CREATE | MPI_MODE_WRONLY,
707 info,
708 &fh);
709 AssertThrowMPI(ierr);
710
711 ierr = MPI_File_set_size(fh, 0); // delete the file contents
712 AssertThrowMPI(ierr);
713 // this barrier is necessary, because otherwise others might already
714 // write while one core is still setting the size to zero.
715 ierr = MPI_Barrier(mpi_communicator);
716 AssertThrowMPI(ierr);
717 ierr = MPI_Info_free(&info);
718 AssertThrowMPI(ierr);
719 // ------------------
720
721 // Write cumulative sizes to file.
722 // Since each processor owns the same information about the data
723 // sizes, it is sufficient to let only the first processor perform
724 // this task.
725 if (myrank == 0)
726 {
728 fh,
729 0,
730 sizes_fixed_cumulative.data(),
731 sizes_fixed_cumulative.size(),
732 MPI_UNSIGNED,
733 MPI_STATUS_IGNORE);
734 AssertThrowMPI(ierr);
735 }
736
737 // Write packed data to file simultaneously.
738 const MPI_Offset size_header =
739 sizes_fixed_cumulative.size() * sizeof(unsigned int);
740
741 // Make sure we do the following computation in 64bit integers to be
742 // able to handle 4GB+ files:
743 const MPI_Offset my_global_file_position =
744 size_header +
745 static_cast<MPI_Offset>(global_first_cell) * bytes_per_cell;
746
747 ierr =
749 my_global_file_position,
750 src_data_fixed.data(),
751 src_data_fixed.size(),
752 MPI_BYTE,
753 MPI_STATUS_IGNORE);
754 AssertThrowMPI(ierr);
755
756 ierr = MPI_File_close(&fh);
757 AssertThrowMPI(ierr);
758 }
759
760
761
762 //
763 // ---------- Variable size data ----------
764 //
765 if (variable_size_data_stored)
766 {
767 const std::string fname_variable =
768 std::string(filename) + "_variable.data";
769
770 MPI_Info info;
771 int ierr = MPI_Info_create(&info);
772 AssertThrowMPI(ierr);
773
774 MPI_File fh;
775 ierr = MPI_File_open(mpi_communicator,
776 fname_variable.c_str(),
777 MPI_MODE_CREATE | MPI_MODE_WRONLY,
778 info,
779 &fh);
780 AssertThrowMPI(ierr);
781
782 ierr = MPI_File_set_size(fh, 0); // delete the file contents
783 AssertThrowMPI(ierr);
784 // this barrier is necessary, because otherwise others might already
785 // write while one core is still setting the size to zero.
786 ierr = MPI_Barrier(mpi_communicator);
787 AssertThrowMPI(ierr);
788 ierr = MPI_Info_free(&info);
789 AssertThrowMPI(ierr);
790
791 // Write sizes of each cell into file simultaneously.
792 {
793 const MPI_Offset my_global_file_position =
794 static_cast<MPI_Offset>(global_first_cell) *
795 sizeof(unsigned int);
796
797 // It is very unlikely that a single process has more than
798 // 2 billion cells, but we might as well check.
799 AssertThrow(src_sizes_variable.size() <
800 static_cast<std::size_t>(
801 std::numeric_limits<int>::max()),
803
805 fh,
806 my_global_file_position,
807 src_sizes_variable.data(),
808 src_sizes_variable.size(),
809 MPI_INT,
810 MPI_STATUS_IGNORE);
811 AssertThrowMPI(ierr);
812 }
813
814 // Gather size of data in bytes we want to store from this
815 // processor and compute the prefix sum. We do this in 64 bit
816 // to avoid overflow for files larger than 4GB:
817 const std::uint64_t size_on_proc = src_data_variable.size();
818 std::uint64_t prefix_sum = 0;
819 ierr = MPI_Exscan(&size_on_proc,
820 &prefix_sum,
821 1,
822 MPI_UINT64_T,
823 MPI_SUM,
824 mpi_communicator);
825 AssertThrowMPI(ierr);
826
827 const MPI_Offset my_global_file_position =
828 static_cast<MPI_Offset>(global_num_cells) * sizeof(unsigned int) +
829 prefix_sum;
830
831 // Write data consecutively into file.
833 fh,
834 my_global_file_position,
835 src_data_variable.data(),
836 src_data_variable.size(),
837 MPI_BYTE,
838 MPI_STATUS_IGNORE);
839 AssertThrowMPI(ierr);
840
841
842 ierr = MPI_File_close(&fh);
843 AssertThrowMPI(ierr);
844 }
845 } // if (mpisize > 1)
846 else
847#endif
848 {
849 (void)global_first_cell;
850 (void)global_num_cells;
851 (void)mpi_communicator;
852
853 //
854 // ---------- Fixed size data ----------
855 //
856 {
857 const std::string fname_fixed = std::string(filename) + "_fixed.data";
858
859 std::ofstream file(fname_fixed, std::ios::binary | std::ios::out);
860
861 // Write header data.
862 file.write(reinterpret_cast<const char *>(
863 sizes_fixed_cumulative.data()),
864 sizes_fixed_cumulative.size() * sizeof(unsigned int));
865
866 // Write packed data.
867 file.write(reinterpret_cast<const char *>(src_data_fixed.data()),
868 src_data_fixed.size() * sizeof(char));
869
870 file.close();
871 }
872
873 //
874 // ---------- Variable size data ----------
875 //
876 if (variable_size_data_stored)
877 {
878 const std::string fname_variable =
879 std::string(filename) + "_variable.data";
880
881 std::ofstream file(fname_variable,
882 std::ios::binary | std::ios::out);
883
884 // Write header data.
885 file.write(reinterpret_cast<const char *>(
886 src_sizes_variable.data()),
887 src_sizes_variable.size() * sizeof(int));
888
889 // Write packed data.
890 file.write(reinterpret_cast<const char *>(src_data_variable.data()),
891 src_data_variable.size() * sizeof(char));
892
893 file.close();
894 }
895 }
896 }
897
898
899 template <int dim, int spacedim>
901 void CellAttachedDataSerializer<dim, spacedim>::load(
902 const unsigned int global_first_cell,
903 const unsigned int global_num_cells,
904 const unsigned int local_num_cells,
905 const std::string &filename,
906 const unsigned int n_attached_deserialize_fixed,
907 const unsigned int n_attached_deserialize_variable,
908 const MPI_Comm &mpi_communicator)
909 {
910 Assert(dest_data_fixed.empty(),
911 ExcMessage("Previously loaded data has not been released yet!"));
912
913 variable_size_data_stored = (n_attached_deserialize_variable > 0);
914
915#ifdef DEAL_II_WITH_MPI
916 // Large fractions of this function have been copied from
917 // DataOutInterface::write_vtu_in_parallel.
918 // TODO: Write general MPIIO interface.
919
920 const int mpisize = Utilities::MPI::n_mpi_processes(mpi_communicator);
921
922 if (mpisize > 1)
923 {
924 //
925 // ---------- Fixed size data ----------
926 //
927 {
928 const std::string fname_fixed = std::string(filename) + "_fixed.data";
929
930 MPI_Info info;
931 int ierr = MPI_Info_create(&info);
932 AssertThrowMPI(ierr);
933
934 MPI_File fh;
935 ierr = MPI_File_open(
936 mpi_communicator, fname_fixed.c_str(), MPI_MODE_RDONLY, info, &fh);
937 AssertThrowMPI(ierr);
938
939 ierr = MPI_Info_free(&info);
940 AssertThrowMPI(ierr);
941
942 // Read cumulative sizes from file.
943 // Since all processors need the same information about the data
944 // sizes, let each of them retrieve it by reading from the same
945 // location in the file.
946 sizes_fixed_cumulative.resize(1 + n_attached_deserialize_fixed +
947 (variable_size_data_stored ? 1 : 0));
949 fh,
950 0,
951 sizes_fixed_cumulative.data(),
952 sizes_fixed_cumulative.size(),
953 MPI_UNSIGNED,
954 MPI_STATUS_IGNORE);
955 AssertThrowMPI(ierr);
956
957 // Allocate sufficient memory.
958 const unsigned int bytes_per_cell = sizes_fixed_cumulative.back();
959 dest_data_fixed.resize(static_cast<size_t>(local_num_cells) *
960 bytes_per_cell);
961
962 // Read packed data from file simultaneously.
963 const MPI_Offset size_header =
964 sizes_fixed_cumulative.size() * sizeof(unsigned int);
965
966 // Make sure we do the following computation in 64bit integers to be
967 // able to handle 4GB+ files:
968 const MPI_Offset my_global_file_position =
969 size_header +
970 static_cast<MPI_Offset>(global_first_cell) * bytes_per_cell;
971
972 ierr =
974 my_global_file_position,
975 dest_data_fixed.data(),
976 dest_data_fixed.size(),
977 MPI_BYTE,
978 MPI_STATUS_IGNORE);
979 AssertThrowMPI(ierr);
980
981
982 ierr = MPI_File_close(&fh);
983 AssertThrowMPI(ierr);
984 }
985
986 //
987 // ---------- Variable size data ----------
988 //
989 if (variable_size_data_stored)
990 {
991 const std::string fname_variable =
992 std::string(filename) + "_variable.data";
993
994 MPI_Info info;
995 int ierr = MPI_Info_create(&info);
996 AssertThrowMPI(ierr);
997
998 MPI_File fh;
999 ierr = MPI_File_open(mpi_communicator,
1000 fname_variable.c_str(),
1001 MPI_MODE_RDONLY,
1002 info,
1003 &fh);
1004 AssertThrowMPI(ierr);
1005
1006 ierr = MPI_Info_free(&info);
1007 AssertThrowMPI(ierr);
1008
1009 // Read sizes of all locally owned cells.
1010 dest_sizes_variable.resize(local_num_cells);
1011
1012 const MPI_Offset my_global_file_position_sizes =
1013 static_cast<MPI_Offset>(global_first_cell) * sizeof(unsigned int);
1014
1016 fh,
1017 my_global_file_position_sizes,
1018 dest_sizes_variable.data(),
1019 dest_sizes_variable.size(),
1020 MPI_INT,
1021 MPI_STATUS_IGNORE);
1022 AssertThrowMPI(ierr);
1023
1024
1025 // Compute my data size in bytes and compute prefix sum. We do this
1026 // in 64 bit to avoid overflow for files larger than 4 GB:
1027 const std::uint64_t size_on_proc =
1028 std::accumulate(dest_sizes_variable.begin(),
1029 dest_sizes_variable.end(),
1030 0ULL);
1031
1032 std::uint64_t prefix_sum = 0;
1033 ierr = MPI_Exscan(&size_on_proc,
1034 &prefix_sum,
1035 1,
1036 MPI_UINT64_T,
1037 MPI_SUM,
1038 mpi_communicator);
1039 AssertThrowMPI(ierr);
1040
1041 const MPI_Offset my_global_file_position =
1042 static_cast<MPI_Offset>(global_num_cells) * sizeof(unsigned int) +
1043 prefix_sum;
1044
1045 dest_data_variable.resize(size_on_proc);
1046
1048 fh,
1049 my_global_file_position,
1050 dest_data_variable.data(),
1051 dest_data_variable.size(),
1052 MPI_BYTE,
1053 MPI_STATUS_IGNORE);
1054 AssertThrowMPI(ierr);
1055
1056 ierr = MPI_File_close(&fh);
1057 AssertThrowMPI(ierr);
1058 }
1059 }
1060 else // if (mpisize > 1)
1061#endif
1062 {
1063 (void)global_first_cell;
1064 (void)global_num_cells;
1065 (void)mpi_communicator;
1066
1067 //
1068 // ---------- Fixed size data ----------
1069 //
1070 {
1071 const std::string fname_fixed = std::string(filename) + "_fixed.data";
1072
1073 std::ifstream file(fname_fixed, std::ios::binary | std::ios::in);
1074 sizes_fixed_cumulative.resize(1 + n_attached_deserialize_fixed +
1075 (variable_size_data_stored ? 1 : 0));
1076
1077 // Read header data.
1078 file.read(reinterpret_cast<char *>(sizes_fixed_cumulative.data()),
1079 sizes_fixed_cumulative.size() * sizeof(unsigned int));
1080
1081 const unsigned int bytes_per_cell = sizes_fixed_cumulative.back();
1082 dest_data_fixed.resize(static_cast<size_t>(local_num_cells) *
1083 bytes_per_cell);
1084
1085 // Read packed data.
1086 file.read(reinterpret_cast<char *>(dest_data_fixed.data()),
1087 dest_data_fixed.size() * sizeof(char));
1088
1089 file.close();
1090 }
1091
1092 //
1093 // ---------- Variable size data ----------
1094 //
1095 if (variable_size_data_stored)
1096 {
1097 const std::string fname_variable =
1098 std::string(filename) + "_variable.data";
1099
1100 std::ifstream file(fname_variable, std::ios::binary | std::ios::in);
1101
1102 // Read header data.
1103 dest_sizes_variable.resize(local_num_cells);
1104 file.read(reinterpret_cast<char *>(dest_sizes_variable.data()),
1105 dest_sizes_variable.size() * sizeof(int));
1106
1107 // Read packed data.
1108 const std::uint64_t size =
1109 std::accumulate(dest_sizes_variable.begin(),
1110 dest_sizes_variable.end(),
1111 0ULL);
1112 dest_data_variable.resize(size);
1113 file.read(reinterpret_cast<char *>(dest_data_variable.data()),
1114 dest_data_variable.size() * sizeof(char));
1115
1116 file.close();
1117 }
1118 }
1119 }
1120
1121
1122 template <int dim, int spacedim>
1124 void CellAttachedDataSerializer<dim, spacedim>::clear()
1125 {
1126 variable_size_data_stored = false;
1127
1128 // free information about data sizes
1129 sizes_fixed_cumulative.clear();
1130 sizes_fixed_cumulative.shrink_to_fit();
1131
1132 // free fixed size transfer data
1133 src_data_fixed.clear();
1134 src_data_fixed.shrink_to_fit();
1135
1136 dest_data_fixed.clear();
1137 dest_data_fixed.shrink_to_fit();
1138
1139 // free variable size transfer data
1140 src_sizes_variable.clear();
1141 src_sizes_variable.shrink_to_fit();
1142
1143 src_data_variable.clear();
1144 src_data_variable.shrink_to_fit();
1145
1146 dest_sizes_variable.clear();
1147 dest_sizes_variable.shrink_to_fit();
1148
1149 dest_data_variable.clear();
1150 dest_data_variable.shrink_to_fit();
1151 }
1152
1153} // namespace internal
1154
1155// anonymous namespace for internal helper functions
1156namespace
1157{
1158 // return whether the given cell is
1159 // patch_level_1, i.e. determine
1160 // whether either all or none of
1161 // its children are further
1162 // refined. this function can only
1163 // be called for non-active cells.
1164 template <int dim, int spacedim>
1165 bool
1166 cell_is_patch_level_1(
1168 {
1169 Assert(cell->is_active() == false, ExcInternalError());
1170
1171 unsigned int n_active_children = 0;
1172 for (unsigned int i = 0; i < cell->n_children(); ++i)
1173 if (cell->child(i)->is_active())
1174 ++n_active_children;
1175
1176 return (n_active_children == 0) ||
1177 (n_active_children == cell->n_children());
1178 }
1179
1180
1181
1182 // return, whether a given @p cell will be
1183 // coarsened, which is the case if all
1184 // children are active and have their coarsen
1185 // flag set. In case only part of the coarsen
1186 // flags are set, remove them.
1187 template <int dim, int spacedim>
1188 bool
1189 cell_will_be_coarsened(
1191 {
1192 // only cells with children should be
1193 // considered for coarsening
1194
1195 if (cell->has_children())
1196 {
1197 unsigned int children_to_coarsen = 0;
1198 const unsigned int n_children = cell->n_children();
1199
1200 for (unsigned int c = 0; c < n_children; ++c)
1201 if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
1202 ++children_to_coarsen;
1203 if (children_to_coarsen == n_children)
1204 return true;
1205 else
1206 for (unsigned int c = 0; c < n_children; ++c)
1207 if (cell->child(c)->is_active())
1208 cell->child(c)->clear_coarsen_flag();
1209 }
1210 // no children, so no coarsening
1211 // possible. however, no children also
1212 // means that this cell will be in the same
1213 // state as if it had children and was
1214 // coarsened. So, what should we return -
1215 // false or true?
1216 // make sure we do not have to do this at
1217 // all...
1218 Assert(cell->has_children(), ExcInternalError());
1219 // ... and then simply return false
1220 return false;
1221 }
1222
1223
1224 // return, whether the face @p face_no of the
1225 // given @p cell will be refined after the
1226 // current refinement step, considering
1227 // refine and coarsen flags and considering
1228 // only those refinemnts that will be caused
1229 // by the neighboring cell.
1230
1231 // this function is used on both active cells
1232 // and cells with children. on cells with
1233 // children it also of interest to know 'how'
1234 // the face will be refined. thus there is an
1235 // additional third argument @p
1236 // expected_face_ref_case returning just
1237 // that. be aware, that this variable will
1238 // only contain useful information if this
1239 // function is called for an active cell.
1240 //
1241 // thus, this is an internal function, users
1242 // should call one of the two alternatives
1243 // following below.
1244 template <int dim, int spacedim>
1245 bool
1246 face_will_be_refined_by_neighbor_internal(
1248 const unsigned int face_no,
1249 RefinementCase<dim - 1> &expected_face_ref_case)
1250 {
1251 // first of all: set the default value for
1252 // expected_face_ref_case, which is no
1253 // refinement at all
1254 expected_face_ref_case = RefinementCase<dim - 1>::no_refinement;
1255
1256 const typename Triangulation<dim, spacedim>::cell_iterator neighbor =
1257 cell->neighbor(face_no);
1258
1259 // If we are at the boundary, there is no
1260 // neighbor which could refine the face
1261 if (neighbor.state() != IteratorState::valid)
1262 return false;
1263
1264 if (neighbor->has_children())
1265 {
1266 // if the neighbor is refined, it may be
1267 // coarsened. if so, then it won't refine
1268 // the face, no matter what else happens
1269 if (cell_will_be_coarsened(neighbor))
1270 return false;
1271 else
1272 // if the neighbor is refined, then it
1273 // is also refined at our current
1274 // face. It will stay so without
1275 // coarsening, so return true in that
1276 // case.
1277 {
1278 expected_face_ref_case = cell->face(face_no)->refinement_case();
1279 return true;
1280 }
1281 }
1282
1283 // now, the neighbor is not refined, but
1284 // perhaps it will be
1285 const RefinementCase<dim> nb_ref_flag = neighbor->refine_flag_set();
1286 if (nb_ref_flag != RefinementCase<dim>::no_refinement)
1287 {
1288 // now we need to know, which of the
1289 // neighbors faces points towards us
1290 const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
1291 // check, whether the cell will be
1292 // refined in a way that refines our
1293 // face
1294 const RefinementCase<dim - 1> face_ref_case =
1296 nb_ref_flag,
1297 neighbor_neighbor,
1298 neighbor->face_orientation(neighbor_neighbor),
1299 neighbor->face_flip(neighbor_neighbor),
1300 neighbor->face_rotation(neighbor_neighbor));
1301 if (face_ref_case != RefinementCase<dim - 1>::no_refinement)
1302 {
1304 neighbor_face = neighbor->face(neighbor_neighbor);
1305 const int this_face_index = cell->face_index(face_no);
1306
1307 // there are still two basic
1308 // possibilities here: the neighbor
1309 // might be coarser or as coarse
1310 // as we are
1311 if (neighbor_face->index() == this_face_index)
1312 // the neighbor is as coarse as
1313 // we are and will be refined at
1314 // the face of consideration, so
1315 // return true
1316 {
1317 expected_face_ref_case = face_ref_case;
1318 return true;
1319 }
1320 else
1321 {
1322 // the neighbor is coarser.
1323 // this is the most complicated
1324 // case. It might be, that the
1325 // neighbor's face will be
1326 // refined, but that we will
1327 // not see this, as we are
1328 // refined in a similar way.
1329
1330 // so, the neighbor's face must
1331 // have children. check, if our
1332 // cell's face is one of these
1333 // (it could also be a
1334 // grand_child)
1335 for (unsigned int c = 0; c < neighbor_face->n_children(); ++c)
1336 if (neighbor_face->child_index(c) == this_face_index)
1337 {
1338 // if the flagged refine
1339 // case of the face is a
1340 // subset or the same as
1341 // the current refine case,
1342 // then the face, as seen
1343 // from our cell, won't be
1344 // refined by the neighbor
1345 if ((neighbor_face->refinement_case() | face_ref_case) ==
1346 neighbor_face->refinement_case())
1347 return false;
1348 else
1349 {
1350 // if we are active, we
1351 // must be an
1352 // anisotropic child
1353 // and the coming
1354 // face_ref_case is
1355 // isotropic. Thus,
1356 // from our cell we
1357 // will see exactly the
1358 // opposite refine case
1359 // that the face has
1360 // now...
1361 Assert(
1362 face_ref_case ==
1365 expected_face_ref_case =
1366 ~neighbor_face->refinement_case();
1367 return true;
1368 }
1369 }
1370
1371 // so, obviously we were not
1372 // one of the children, but a
1373 // grandchild. This is only
1374 // possible in 3d.
1375 Assert(dim == 3, ExcInternalError());
1376 // In that case, however, no
1377 // matter what the neighbor
1378 // does, it won't be finer
1379 // after the next refinement
1380 // step.
1381 return false;
1382 }
1383 } // if face will be refined
1384 } // if neighbor is flagged for refinement
1385
1386 // no cases left, so the neighbor will not
1387 // refine the face
1388 return false;
1389 }
1390
1391 // version of above function for both active
1392 // and non-active cells
1393 template <int dim, int spacedim>
1394 bool
1395 face_will_be_refined_by_neighbor(
1397 const unsigned int face_no)
1398 {
1399 RefinementCase<dim - 1> dummy = RefinementCase<dim - 1>::no_refinement;
1400 return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
1401 }
1402
1403 // version of above function for active cells
1404 // only. Additionally returning the refine
1405 // case (to come) of the face under
1406 // consideration
1407 template <int dim, int spacedim>
1408 bool
1409 face_will_be_refined_by_neighbor(
1411 const unsigned int face_no,
1412 RefinementCase<dim - 1> &expected_face_ref_case)
1413 {
1414 return face_will_be_refined_by_neighbor_internal(cell,
1415 face_no,
1416 expected_face_ref_case);
1417 }
1418
1419
1420
1421 template <int dim, int spacedim>
1422 bool
1423 satisfies_level1_at_vertex_rule(
1425 {
1426 std::vector<unsigned int> min_adjacent_cell_level(
1427 triangulation.n_vertices(), triangulation.n_levels());
1428 std::vector<unsigned int> max_adjacent_cell_level(
1429 triangulation.n_vertices(), 0);
1430
1431 for (const auto &cell : triangulation.active_cell_iterators())
1432 for (const unsigned int v : cell->vertex_indices())
1433 {
1434 min_adjacent_cell_level[cell->vertex_index(v)] =
1435 std::min<unsigned int>(
1436 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
1437 max_adjacent_cell_level[cell->vertex_index(v)] =
1438 std::max<unsigned int>(
1439 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
1440 }
1441
1442 for (unsigned int k = 0; k < triangulation.n_vertices(); ++k)
1443 if (triangulation.vertex_used(k))
1444 if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
1445 return false;
1446 return true;
1447 }
1448
1449
1450
1468 template <int dim, int spacedim>
1469 unsigned int
1470 middle_vertex_index(
1472 {
1473 if (line->has_children())
1474 return line->child(0)->vertex_index(1);
1476 }
1477
1478
1479 template <int dim, int spacedim>
1480 unsigned int
1481 middle_vertex_index(
1483 {
1484 switch (static_cast<unsigned char>(quad->refinement_case()))
1485 {
1487 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
1488 break;
1490 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
1491 break;
1493 return quad->child(0)->vertex_index(3);
1494 break;
1495 default:
1496 break;
1497 }
1499 }
1500
1501
1502 template <int dim, int spacedim>
1503 unsigned int
1504 middle_vertex_index(
1506 {
1507 switch (static_cast<unsigned char>(hex->refinement_case()))
1508 {
1510 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
1511 break;
1513 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
1514 break;
1516 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
1517 break;
1519 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
1520 break;
1522 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
1523 break;
1525 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
1526 break;
1528 return hex->child(0)->vertex_index(7);
1529 break;
1530 default:
1531 break;
1532 }
1534 }
1535
1536
1549 template <class TRIANGULATION>
1550 inline typename TRIANGULATION::DistortedCellList
1551 collect_distorted_coarse_cells(const TRIANGULATION &)
1552 {
1553 return typename TRIANGULATION::DistortedCellList();
1554 }
1555
1556
1557
1566 template <int dim>
1568 collect_distorted_coarse_cells(const Triangulation<dim, dim> &triangulation)
1569 {
1570 typename Triangulation<dim, dim>::DistortedCellList distorted_cells;
1571 for (const auto &cell : triangulation.cell_iterators_on_level(0))
1572 {
1574 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1575 vertices[i] = cell->vertex(i);
1576
1579
1580 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1581 if (determinants[i] <=
1582 1e-9 * Utilities::fixed_power<dim>(cell->diameter()))
1583 {
1584 distorted_cells.distorted_cells.push_back(cell);
1585 break;
1586 }
1587 }
1588
1589 return distorted_cells;
1590 }
1591
1592
1599 template <int dim>
1600 bool
1601 has_distorted_children(
1602 const typename Triangulation<dim, dim>::cell_iterator &cell)
1603 {
1604 Assert(cell->has_children(), ExcInternalError());
1605
1606 for (unsigned int c = 0; c < cell->n_children(); ++c)
1607 {
1609 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1610 vertices[i] = cell->child(c)->vertex(i);
1611
1614
1615 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1616 if (determinants[i] <=
1617 1e-9 * Utilities::fixed_power<dim>(cell->child(c)->diameter()))
1618 return true;
1619 }
1620
1621 return false;
1622 }
1623
1624
1632 template <int dim, int spacedim>
1633 bool
1634 has_distorted_children(
1636 {
1637 return false;
1638 }
1639
1640
1641 template <int dim, int spacedim>
1642 void
1643 update_periodic_face_map_recursively(
1644 const typename Triangulation<dim, spacedim>::cell_iterator &cell_1,
1645 const typename Triangulation<dim, spacedim>::cell_iterator &cell_2,
1646 unsigned int n_face_1,
1647 unsigned int n_face_2,
1648 const unsigned char orientation,
1649 typename std::map<
1651 unsigned int>,
1652 std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
1653 unsigned int>,
1654 unsigned char>> &periodic_face_map)
1655 {
1656 using FaceIterator = typename Triangulation<dim, spacedim>::face_iterator;
1657 const FaceIterator face_1 = cell_1->face(n_face_1);
1658 const FaceIterator face_2 = cell_2->face(n_face_2);
1659
1660 const unsigned char inverse_orientation =
1661 face_1->reference_cell().get_inverse_combined_orientation(orientation);
1662
1663#ifdef DEBUG
1664 const auto [face_orientation, face_rotation, face_flip] =
1666
1667 Assert((dim != 1) || (face_orientation == true && face_flip == false &&
1668 face_rotation == false),
1669 ExcMessage("The supplied orientation "
1670 "(face_orientation, face_flip, face_rotation) "
1671 "is invalid for 1d"));
1672
1673 Assert((dim != 2) || (face_flip == false && face_rotation == false),
1674 ExcMessage("The supplied orientation "
1675 "(face_orientation, face_flip, face_rotation) "
1676 "is invalid for 2d"));
1677#endif
1678
1679 Assert(face_1 != face_2, ExcMessage("face_1 and face_2 are equal!"));
1680
1681 Assert(face_1->at_boundary() && face_2->at_boundary(),
1682 ExcMessage("Periodic faces must be on the boundary"));
1683
1684 // Check if the requirement that each edge can only have at most one hanging
1685 // node, and as a consequence neighboring cells can differ by at most
1686 // one refinement level is enforced. In 1d, there are no hanging nodes and
1687 // so neighboring cells can differ by more than one refinement level.
1688 Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2,
1690
1691 // insert periodic face pair for both cells
1692 using CellFace =
1693 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
1694 unsigned int>;
1695 const CellFace cell_face_1(cell_1, n_face_1);
1696 const CellFace cell_face_2(cell_2, n_face_2);
1697 const std::pair<CellFace, unsigned char> cell_face_orientation_2(
1698 cell_face_2, orientation);
1699
1700 const std::pair<CellFace, std::pair<CellFace, unsigned char>>
1701 periodic_faces(cell_face_1, cell_face_orientation_2);
1702
1703 // Only one periodic neighbor is allowed
1704 Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError());
1705 periodic_face_map.insert(periodic_faces);
1706
1707 if (dim == 1)
1708 {
1709 if (cell_1->has_children())
1710 {
1711 if (cell_2->has_children())
1712 {
1713 update_periodic_face_map_recursively<dim, spacedim>(
1714 cell_1->child(n_face_1),
1715 cell_2->child(n_face_2),
1716 n_face_1,
1717 n_face_2,
1718 orientation,
1719 periodic_face_map);
1720 }
1721 else // only face_1 has children
1722 {
1723 update_periodic_face_map_recursively<dim, spacedim>(
1724 cell_1->child(n_face_1),
1725 cell_2,
1726 n_face_1,
1727 n_face_2,
1728 orientation,
1729 periodic_face_map);
1730 }
1731 }
1732 }
1733 else // dim == 2 || dim == 3
1734 {
1735 if (cell_1->has_children())
1736 {
1737 if (cell_2->has_children())
1738 {
1739 // In the case that both faces have children, we loop over all
1740 // children and apply update_periodic_face_map_recursively
1741 // recursively:
1742
1743 Assert(face_1->n_children() ==
1745 face_2->n_children() ==
1748
1749 const auto reference_cell = cell_1->reference_cell();
1750
1751 for (unsigned int i = 0;
1752 i < GeometryInfo<dim>::max_children_per_face;
1753 ++i)
1754 {
1755 // Lookup the index for the second face
1756 const unsigned int j =
1757 reference_cell.standard_to_real_face_vertex(
1758 i, n_face_1, inverse_orientation);
1759
1760 // find subcell ids that belong to the subface indices
1761 unsigned int child_cell_1 =
1763 cell_1->refinement_case(),
1764 n_face_1,
1765 i,
1766 cell_1->face_orientation(n_face_1),
1767 cell_1->face_flip(n_face_1),
1768 cell_1->face_rotation(n_face_1),
1769 face_1->refinement_case());
1770 unsigned int child_cell_2 =
1772 cell_2->refinement_case(),
1773 n_face_2,
1774 j,
1775 cell_2->face_orientation(n_face_2),
1776 cell_2->face_flip(n_face_2),
1777 cell_2->face_rotation(n_face_2),
1778 face_2->refinement_case());
1779
1780 Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
1781 face_1->child(i),
1783 Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
1784 face_2->child(j),
1786
1787 // precondition: subcell has the same orientation as cell
1788 // (so that the face numbers coincide) recursive call
1789 update_periodic_face_map_recursively<dim, spacedim>(
1790 cell_1->child(child_cell_1),
1791 cell_2->child(child_cell_2),
1792 n_face_1,
1793 n_face_2,
1794 orientation,
1795 periodic_face_map);
1796 }
1797 }
1798 else // only face_1 has children
1799 {
1800 for (unsigned int i = 0;
1801 i < GeometryInfo<dim>::max_children_per_face;
1802 ++i)
1803 {
1804 // find subcell ids that belong to the subface indices
1805 unsigned int child_cell_1 =
1807 cell_1->refinement_case(),
1808 n_face_1,
1809 i,
1810 cell_1->face_orientation(n_face_1),
1811 cell_1->face_flip(n_face_1),
1812 cell_1->face_rotation(n_face_1),
1813 face_1->refinement_case());
1814
1815 // recursive call
1816 update_periodic_face_map_recursively<dim, spacedim>(
1817 cell_1->child(child_cell_1),
1818 cell_2,
1819 n_face_1,
1820 n_face_2,
1821 orientation,
1822 periodic_face_map);
1823 }
1824 }
1825 }
1826 }
1827 }
1828
1829
1830} // end of anonymous namespace
1831
1832
1833namespace internal
1834{
1835 namespace TriangulationImplementation
1836 {
1837 // make sure that if in the following we
1838 // write Triangulation<dim,spacedim>
1839 // we mean the *class*
1840 // ::Triangulation, not the
1841 // enclosing namespace
1842 // internal::TriangulationImplementation
1843 using ::Triangulation;
1844
1850 int,
1851 << "Something went wrong upon construction of cell "
1852 << arg1);
1863 int,
1864 << "Cell " << arg1
1865 << " has negative measure. This typically "
1866 << "indicates some distortion in the cell, or a mistakenly "
1867 << "swapped pair of vertices in the input to "
1868 << "Triangulation::create_triangulation().");
1877 int,
1878 int,
1879 int,
1880 << "Error while creating cell " << arg1
1881 << ": the vertex index " << arg2 << " must be between 0 and "
1882 << arg3 << '.');
1889 int,
1890 int,
1892 << "The input data for creating a triangulation contained "
1893 << "information about a line with indices " << arg1 << " and " << arg2
1894 << " that is described to have boundary indicator "
1895 << static_cast<int>(arg3)
1896 << ". However, this is an internal line not located on the "
1897 << "boundary. You cannot assign a boundary indicator to it." << std::endl
1898 << std::endl
1899 << "If this happened at a place where you call "
1900 << "Triangulation::create_triangulation() yourself, you need "
1901 << "to check the SubCellData object you pass to this function."
1902 << std::endl
1903 << std::endl
1904 << "If this happened in a place where you are reading a mesh "
1905 << "from a file, then you need to investigate why such a line "
1906 << "ended up in the input file. A typical case is a geometry "
1907 << "that consisted of multiple parts and for which the mesh "
1908 << "generator program assumes that the interface between "
1909 << "two parts is a boundary when that isn't supposed to be "
1910 << "the case, or where the mesh generator simply assigns "
1911 << "'geometry indicators' to lines at the perimeter of "
1912 << "a part that are not supposed to be interpreted as "
1913 << "'boundary indicators'.");
1920 int,
1921 int,
1922 int,
1923 int,
1925 << "The input data for creating a triangulation contained "
1926 << "information about a quad with indices " << arg1 << ", " << arg2
1927 << ", " << arg3 << ", and " << arg4
1928 << " that is described to have boundary indicator "
1929 << static_cast<int>(arg5)
1930 << ". However, this is an internal quad not located on the "
1931 << "boundary. You cannot assign a boundary indicator to it." << std::endl
1932 << std::endl
1933 << "If this happened at a place where you call "
1934 << "Triangulation::create_triangulation() yourself, you need "
1935 << "to check the SubCellData object you pass to this function."
1936 << std::endl
1937 << std::endl
1938 << "If this happened in a place where you are reading a mesh "
1939 << "from a file, then you need to investigate why such a quad "
1940 << "ended up in the input file. A typical case is a geometry "
1941 << "that consisted of multiple parts and for which the mesh "
1942 << "generator program assumes that the interface between "
1943 << "two parts is a boundary when that isn't supposed to be "
1944 << "the case, or where the mesh generator simply assigns "
1945 << "'geometry indicators' to quads at the surface of "
1946 << "a part that are not supposed to be interpreted as "
1947 << "'boundary indicators'.");
1954 int,
1955 int,
1956 << "In SubCellData the line info of the line with vertex indices " << arg1
1957 << " and " << arg2 << " appears more than once. "
1958 << "This is not allowed.");
1965 int,
1966 int,
1967 std::string,
1968 << "In SubCellData the line info of the line with vertex indices " << arg1
1969 << " and " << arg2 << " appears multiple times with different (valid) "
1970 << arg3 << ". This is not allowed.");
1977 int,
1978 int,
1979 int,
1980 int,
1981 std::string,
1982 << "In SubCellData the quad info of the quad with line indices " << arg1
1983 << ", " << arg2 << ", " << arg3 << " and " << arg4
1984 << " appears multiple times with different (valid) " << arg5
1985 << ". This is not allowed.");
1986
1987 /*
1988 * Reserve space for TriaFaces. Details:
1989 *
1990 * Reserve space for line_orientations.
1991 *
1992 * @note Used only for dim=3.
1993 */
1994 void
1996 const unsigned int new_quads_in_pairs,
1997 const unsigned int new_quads_single)
1998 {
1999 AssertDimension(tria_faces.dim, 3);
2000
2001 Assert(new_quads_in_pairs % 2 == 0, ExcInternalError());
2002
2003 unsigned int next_free_single = 0;
2004 unsigned int next_free_pair = 0;
2005
2006 // count the number of objects, of unused single objects and of
2007 // unused pairs of objects
2008 unsigned int n_quads = 0;
2009 unsigned int n_unused_pairs = 0;
2010 unsigned int n_unused_singles = 0;
2011 for (unsigned int i = 0; i < tria_faces.quads.used.size(); ++i)
2012 {
2013 if (tria_faces.quads.used[i])
2014 ++n_quads;
2015 else if (i + 1 < tria_faces.quads.used.size())
2016 {
2017 if (tria_faces.quads.used[i + 1])
2018 {
2019 ++n_unused_singles;
2020 if (next_free_single == 0)
2021 next_free_single = i;
2022 }
2023 else
2024 {
2025 ++n_unused_pairs;
2026 if (next_free_pair == 0)
2027 next_free_pair = i;
2028 ++i;
2029 }
2030 }
2031 else
2032 ++n_unused_singles;
2033 }
2034 Assert(n_quads + 2 * n_unused_pairs + n_unused_singles ==
2035 tria_faces.quads.used.size(),
2037 (void)n_quads;
2038
2039 // how many single quads are needed in addition to n_unused_quads?
2040 const int additional_single_quads = new_quads_single - n_unused_singles;
2041
2042 unsigned int new_size =
2043 tria_faces.quads.used.size() + new_quads_in_pairs - 2 * n_unused_pairs;
2044 if (additional_single_quads > 0)
2045 new_size += additional_single_quads;
2046
2047 // see above...
2048 if (new_size > tria_faces.quads.n_objects())
2049 {
2050 // reserve the field of the derived class
2051 tria_faces.quads_line_orientations.resize(
2052 new_size * GeometryInfo<3>::lines_per_face, true);
2053
2054 auto &q_is_q = tria_faces.quad_is_quadrilateral;
2055 q_is_q.reserve(new_size);
2056 q_is_q.insert(q_is_q.end(), new_size - q_is_q.size(), true);
2057 }
2058 }
2059
2060
2061
2075 void
2077 const unsigned int total_cells,
2078 const unsigned int dimension,
2079 const unsigned int space_dimension,
2080 const bool tetraheder_in_mesh = false)
2081 {
2082 // we need space for total_cells cells. Maybe we have more already
2083 // with those cells which are unused, so only allocate new space if
2084 // needed.
2085 //
2086 // note that all arrays should have equal sizes (checked by
2087 // @p{monitor_memory}
2088 if (total_cells > tria_level.refine_flags.size())
2089 {
2090 tria_level.refine_flags.reserve(total_cells);
2091 tria_level.refine_flags.insert(tria_level.refine_flags.end(),
2092 total_cells -
2093 tria_level.refine_flags.size(),
2094 /*RefinementCase::no_refinement=*/0);
2095
2096 if (tetraheder_in_mesh)
2097 {
2098 tria_level.refine_choice.reserve(total_cells);
2099 tria_level.refine_choice.insert(
2100 tria_level.refine_choice.end(),
2101 total_cells - tria_level.refine_choice.size(),
2102 static_cast<char>(
2104 }
2105
2106 tria_level.coarsen_flags.reserve(total_cells);
2107 tria_level.coarsen_flags.insert(tria_level.coarsen_flags.end(),
2108 total_cells -
2109 tria_level.coarsen_flags.size(),
2110 false);
2111
2112 tria_level.active_cell_indices.reserve(total_cells);
2113 tria_level.active_cell_indices.insert(
2114 tria_level.active_cell_indices.end(),
2115 total_cells - tria_level.active_cell_indices.size(),
2117
2118 tria_level.subdomain_ids.reserve(total_cells);
2119 tria_level.subdomain_ids.insert(tria_level.subdomain_ids.end(),
2120 total_cells -
2121 tria_level.subdomain_ids.size(),
2122 0);
2123
2124 tria_level.level_subdomain_ids.reserve(total_cells);
2125 tria_level.level_subdomain_ids.insert(
2126 tria_level.level_subdomain_ids.end(),
2127 total_cells - tria_level.level_subdomain_ids.size(),
2128 0);
2129
2130 tria_level.global_active_cell_indices.reserve(total_cells);
2131 tria_level.global_active_cell_indices.insert(
2132 tria_level.global_active_cell_indices.end(),
2133 total_cells - tria_level.global_active_cell_indices.size(),
2135
2136 tria_level.global_level_cell_indices.reserve(total_cells);
2137 tria_level.global_level_cell_indices.insert(
2138 tria_level.global_level_cell_indices.end(),
2139 total_cells - tria_level.global_level_cell_indices.size(),
2141
2142 if (dimension == space_dimension - 1)
2143 {
2144 tria_level.direction_flags.reserve(total_cells);
2145 tria_level.direction_flags.insert(
2146 tria_level.direction_flags.end(),
2147 total_cells - tria_level.direction_flags.size(),
2148 true);
2149 }
2150 else
2151 tria_level.direction_flags.clear();
2152
2153 tria_level.parents.reserve((total_cells + 1) / 2);
2154 tria_level.parents.insert(tria_level.parents.end(),
2155 (total_cells + 1) / 2 -
2156 tria_level.parents.size(),
2157 -1);
2158
2159 tria_level.neighbors.reserve(total_cells * (2 * dimension));
2160 tria_level.neighbors.insert(tria_level.neighbors.end(),
2161 total_cells * (2 * dimension) -
2162 tria_level.neighbors.size(),
2163 std::make_pair(-1, -1));
2164
2165 if (tria_level.dim == 2 || tria_level.dim == 3)
2166 {
2167 const unsigned int max_faces_per_cell = 2 * dimension;
2168 tria_level.face_orientations.resize(total_cells *
2169 max_faces_per_cell);
2170
2171 tria_level.reference_cell.reserve(total_cells);
2172 tria_level.reference_cell.insert(
2173 tria_level.reference_cell.end(),
2174 total_cells - tria_level.reference_cell.size(),
2175 tria_level.dim == 2 ? ReferenceCells::Quadrilateral :
2177 }
2178 }
2179 }
2180
2181
2182
2187 int,
2188 int,
2189 << "The containers have sizes " << arg1 << " and " << arg2
2190 << ", which is not as expected.");
2191
2197 void
2198 monitor_memory(const TriaLevel &tria_level,
2199 const unsigned int true_dimension)
2200 {
2201 (void)tria_level;
2202 (void)true_dimension;
2203 Assert(2 * true_dimension * tria_level.refine_flags.size() ==
2204 tria_level.neighbors.size(),
2205 ExcMemoryInexact(tria_level.refine_flags.size(),
2206 tria_level.neighbors.size()));
2207 Assert(2 * true_dimension * tria_level.coarsen_flags.size() ==
2208 tria_level.neighbors.size(),
2209 ExcMemoryInexact(tria_level.coarsen_flags.size(),
2210 tria_level.neighbors.size()));
2211 }
2212
2213
2214
2227 void
2229 const unsigned int new_objects_in_pairs,
2230 const unsigned int new_objects_single = 0)
2231 {
2232 if (tria_objects.structdim <= 2)
2233 {
2234 Assert(new_objects_in_pairs % 2 == 0, ExcInternalError());
2235
2236 tria_objects.next_free_single = 0;
2237 tria_objects.next_free_pair = 0;
2238 tria_objects.reverse_order_next_free_single = false;
2239
2240 // count the number of objects, of unused single objects and of
2241 // unused pairs of objects
2242 unsigned int n_objects = 0;
2243 unsigned int n_unused_pairs = 0;
2244 unsigned int n_unused_singles = 0;
2245 for (unsigned int i = 0; i < tria_objects.used.size(); ++i)
2246 {
2247 if (tria_objects.used[i])
2248 ++n_objects;
2249 else if (i + 1 < tria_objects.used.size())
2250 {
2251 if (tria_objects.used[i + 1])
2252 {
2253 ++n_unused_singles;
2254 if (tria_objects.next_free_single == 0)
2255 tria_objects.next_free_single = i;
2256 }
2257 else
2258 {
2259 ++n_unused_pairs;
2260 if (tria_objects.next_free_pair == 0)
2261 tria_objects.next_free_pair = i;
2262 ++i;
2263 }
2264 }
2265 else
2266 ++n_unused_singles;
2267 }
2268 Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
2269 tria_objects.used.size(),
2271 (void)n_objects;
2272
2273 // how many single objects are needed in addition to
2274 // n_unused_objects?
2275 const int additional_single_objects =
2276 new_objects_single - n_unused_singles;
2277
2278 unsigned int new_size = tria_objects.used.size() +
2279 new_objects_in_pairs - 2 * n_unused_pairs;
2280 if (additional_single_objects > 0)
2281 new_size += additional_single_objects;
2282
2283 // only allocate space if necessary
2284 if (new_size > tria_objects.n_objects())
2285 {
2286 const unsigned int max_faces_per_cell =
2287 2 * tria_objects.structdim;
2288 const unsigned int max_children_per_cell =
2289 1 << tria_objects.structdim;
2290
2291 tria_objects.cells.reserve(new_size * max_faces_per_cell);
2292 tria_objects.cells.insert(tria_objects.cells.end(),
2293 (new_size - tria_objects.n_objects()) *
2294 max_faces_per_cell,
2295 -1);
2296
2297 tria_objects.used.reserve(new_size);
2298 tria_objects.used.insert(tria_objects.used.end(),
2299 new_size - tria_objects.used.size(),
2300 false);
2301
2302 tria_objects.user_flags.reserve(new_size);
2303 tria_objects.user_flags.insert(tria_objects.user_flags.end(),
2304 new_size -
2305 tria_objects.user_flags.size(),
2306 false);
2307
2308 const unsigned int factor = max_children_per_cell / 2;
2309 tria_objects.children.reserve(factor * new_size);
2310 tria_objects.children.insert(tria_objects.children.end(),
2311 factor * new_size -
2312 tria_objects.children.size(),
2313 -1);
2314
2315 if (tria_objects.structdim > 1)
2316 {
2317 tria_objects.refinement_cases.reserve(new_size);
2318 tria_objects.refinement_cases.insert(
2319 tria_objects.refinement_cases.end(),
2320 new_size - tria_objects.refinement_cases.size(),
2321 /*RefinementCase::no_refinement=*/0);
2322 }
2323
2324 // first reserve, then resize. Otherwise the std library can
2325 // decide to allocate more entries.
2326 tria_objects.boundary_or_material_id.reserve(new_size);
2327 tria_objects.boundary_or_material_id.resize(new_size);
2328
2329 tria_objects.user_data.reserve(new_size);
2330 tria_objects.user_data.resize(new_size);
2331
2332 tria_objects.manifold_id.reserve(new_size);
2333 tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
2334 new_size -
2335 tria_objects.manifold_id.size(),
2337 }
2338
2339 if (n_unused_singles == 0)
2340 {
2341 tria_objects.next_free_single = new_size - 1;
2342 tria_objects.reverse_order_next_free_single = true;
2343 }
2344 }
2345 else
2346 {
2347 const unsigned int new_hexes = new_objects_in_pairs;
2348
2349 const unsigned int new_size =
2350 new_hexes + std::count(tria_objects.used.begin(),
2351 tria_objects.used.end(),
2352 true);
2353
2354 // see above...
2355 if (new_size > tria_objects.n_objects())
2356 {
2357 const unsigned int max_faces_per_cell =
2358 2 * tria_objects.structdim;
2359
2360 tria_objects.cells.reserve(new_size * max_faces_per_cell);
2361 tria_objects.cells.insert(tria_objects.cells.end(),
2362 (new_size - tria_objects.n_objects()) *
2363 max_faces_per_cell,
2364 -1);
2365
2366 tria_objects.used.reserve(new_size);
2367 tria_objects.used.insert(tria_objects.used.end(),
2368 new_size - tria_objects.used.size(),
2369 false);
2370
2371 tria_objects.user_flags.reserve(new_size);
2372 tria_objects.user_flags.insert(tria_objects.user_flags.end(),
2373 new_size -
2374 tria_objects.user_flags.size(),
2375 false);
2376
2377 tria_objects.children.reserve(4 * new_size);
2378 tria_objects.children.insert(tria_objects.children.end(),
2379 4 * new_size -
2380 tria_objects.children.size(),
2381 -1);
2382
2383 // for the following fields, we know exactly how many elements
2384 // we need, so first reserve then resize (resize itself, at least
2385 // with some compiler libraries, appears to round up the size it
2386 // actually reserves)
2387 tria_objects.boundary_or_material_id.reserve(new_size);
2388 tria_objects.boundary_or_material_id.resize(new_size);
2389
2390 tria_objects.manifold_id.reserve(new_size);
2391 tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
2392 new_size -
2393 tria_objects.manifold_id.size(),
2395
2396 tria_objects.user_data.reserve(new_size);
2397 tria_objects.user_data.resize(new_size);
2398
2399 tria_objects.refinement_cases.reserve(new_size);
2400 tria_objects.refinement_cases.insert(
2401 tria_objects.refinement_cases.end(),
2402 new_size - tria_objects.refinement_cases.size(),
2403 /*RefinementCase::no_refinement=*/0);
2404 }
2405 tria_objects.next_free_single = tria_objects.next_free_pair = 0;
2406 }
2407 }
2408
2409
2410
2416 void
2417 monitor_memory(const TriaObjects &tria_object, const unsigned int)
2418 {
2419 Assert(tria_object.n_objects() == tria_object.used.size(),
2420 ExcMemoryInexact(tria_object.n_objects(),
2421 tria_object.used.size()));
2422 Assert(tria_object.n_objects() == tria_object.user_flags.size(),
2423 ExcMemoryInexact(tria_object.n_objects(),
2424 tria_object.user_flags.size()));
2425 Assert(tria_object.n_objects() ==
2426 tria_object.boundary_or_material_id.size(),
2427 ExcMemoryInexact(tria_object.n_objects(),
2428 tria_object.boundary_or_material_id.size()));
2429 Assert(tria_object.n_objects() == tria_object.manifold_id.size(),
2430 ExcMemoryInexact(tria_object.n_objects(),
2431 tria_object.manifold_id.size()));
2432 Assert(tria_object.n_objects() == tria_object.user_data.size(),
2433 ExcMemoryInexact(tria_object.n_objects(),
2434 tria_object.user_data.size()));
2435
2436 if (tria_object.structdim == 1)
2437 {
2438 Assert(1 * tria_object.n_objects() == tria_object.children.size(),
2439 ExcMemoryInexact(tria_object.n_objects(),
2440 tria_object.children.size()));
2441 }
2442 else if (tria_object.structdim == 2)
2443 {
2444 Assert(2 * tria_object.n_objects() == tria_object.children.size(),
2445 ExcMemoryInexact(tria_object.n_objects(),
2446 tria_object.children.size()));
2447 }
2448 else if (tria_object.structdim == 3)
2449 {
2450 Assert(4 * tria_object.n_objects() == tria_object.children.size(),
2451 ExcMemoryInexact(tria_object.n_objects(),
2452 tria_object.children.size()));
2453 }
2454 }
2455
2456
2457
2462 template <int dim, int spacedim>
2464 {
2465 public:
2469 virtual ~Policy() = default;
2470
2474 virtual void
2476
2480 virtual void
2484 std::vector<unsigned int> &line_cell_count,
2485 std::vector<unsigned int> &quad_cell_count) = 0;
2486
2492 const bool check_for_distorted_cells) = 0;
2493
2497 virtual void
2500
2504 virtual void
2507
2511 virtual bool
2513 const typename Triangulation<dim, spacedim>::cell_iterator &cell) = 0;
2514
2521 virtual std::unique_ptr<Policy<dim, spacedim>>
2522 clone() = 0;
2523 };
2524
2525
2526
2532 template <int dim, int spacedim, typename T>
2533 class PolicyWrapper : public Policy<dim, spacedim>
2534 {
2535 public:
2536 void
2538 {
2539 T::update_neighbors(tria);
2540 }
2541
2542 void
2546 std::vector<unsigned int> &line_cell_count,
2547 std::vector<unsigned int> &quad_cell_count) override
2548 {
2549 T::delete_children(tria, cell, line_cell_count, quad_cell_count);
2550 }
2551
2554 const bool check_for_distorted_cells) override
2555 {
2556 return T::execute_refinement(triangulation, check_for_distorted_cells);
2557 }
2558
2559 void
2562 {
2563 T::prevent_distorted_boundary_cells(triangulation);
2564 }
2565
2566 void
2569 {
2570 T::prepare_refinement_dim_dependent(triangulation);
2571 }
2572
2573 bool
2576 override
2577 {
2578 return T::template coarsening_allowed<dim, spacedim>(cell);
2579 }
2580
2581 std::unique_ptr<Policy<dim, spacedim>>
2582 clone() override
2583 {
2584 return std::make_unique<PolicyWrapper<dim, spacedim, T>>();
2585 }
2586 };
2587
2588
2589
2686 {
2698 template <int dim, int spacedim>
2699 static void
2702 const unsigned int level_objects,
2704 {
2705 using line_iterator =
2707
2708 number_cache.n_levels = 0;
2709 if (level_objects > 0)
2710 // find the last level on which there are used cells
2711 for (unsigned int level = 0; level < level_objects; ++level)
2712 if (triangulation.begin(level) != triangulation.end(level))
2713 number_cache.n_levels = level + 1;
2714
2715 // no cells at all?
2716 Assert(number_cache.n_levels > 0, ExcInternalError());
2717
2718 //---------------------------------
2719 // update the number of lines on the different levels in the
2720 // cache
2721 number_cache.n_lines = 0;
2722 number_cache.n_active_lines = 0;
2723
2724 // for 1d, lines have levels so take count the objects per
2725 // level and globally
2726 if (dim == 1)
2727 {
2728 number_cache.n_lines_level.resize(number_cache.n_levels);
2729 number_cache.n_active_lines_level.resize(number_cache.n_levels);
2730
2731 for (unsigned int level = 0; level < number_cache.n_levels; ++level)
2732 {
2733 // count lines on this level
2734 number_cache.n_lines_level[level] = 0;
2735 number_cache.n_active_lines_level[level] = 0;
2736
2737 line_iterator line = triangulation.begin_line(level),
2738 endc =
2739 (level == number_cache.n_levels - 1 ?
2740 line_iterator(triangulation.end_line()) :
2741 triangulation.begin_line(level + 1));
2742 for (; line != endc; ++line)
2743 {
2744 ++number_cache.n_lines_level[level];
2745 if (line->has_children() == false)
2746 ++number_cache.n_active_lines_level[level];
2747 }
2748
2749 // update total number of lines
2750 number_cache.n_lines += number_cache.n_lines_level[level];
2751 number_cache.n_active_lines +=
2752 number_cache.n_active_lines_level[level];
2753 }
2754 }
2755 else
2756 {
2757 // for dim>1, there are no levels for lines
2758 number_cache.n_lines_level.clear();
2759 number_cache.n_active_lines_level.clear();
2760
2761 line_iterator line = triangulation.begin_line(),
2762 endc = triangulation.end_line();
2763 for (; line != endc; ++line)
2764 {
2765 ++number_cache.n_lines;
2766 if (line->has_children() == false)
2767 ++number_cache.n_active_lines;
2768 }
2769 }
2770 }
2771
2786 template <int dim, int spacedim>
2787 static void
2790 const unsigned int level_objects,
2792 {
2793 // update lines and n_levels in number_cache. since we don't
2794 // access any of these numbers, we can do this in the
2795 // background
2797 static_cast<
2798 void (*)(const Triangulation<dim, spacedim> &,
2799 const unsigned int,
2801 &compute_number_cache_dim<dim, spacedim>),
2803 level_objects,
2805 number_cache));
2806
2807 using quad_iterator =
2809
2810 //---------------------------------
2811 // update the number of quads on the different levels in the
2812 // cache
2813 number_cache.n_quads = 0;
2814 number_cache.n_active_quads = 0;
2815
2816 // for 2d, quads have levels so take count the objects per
2817 // level and globally
2818 if (dim == 2)
2819 {
2820 // count the number of levels; the function we called above
2821 // on a separate Task for lines also does this and puts it into
2822 // number_cache.n_levels, but this datum may not yet be
2823 // available as we call the function on a separate task
2824 unsigned int n_levels = 0;
2825 if (level_objects > 0)
2826 // find the last level on which there are used cells
2827 for (unsigned int level = 0; level < level_objects; ++level)
2828 if (triangulation.begin(level) != triangulation.end(level))
2829 n_levels = level + 1;
2830
2831 number_cache.n_quads_level.resize(n_levels);
2832 number_cache.n_active_quads_level.resize(n_levels);
2833
2834 for (unsigned int level = 0; level < n_levels; ++level)
2835 {
2836 // count quads on this level
2837 number_cache.n_quads_level[level] = 0;
2838 number_cache.n_active_quads_level[level] = 0;
2839
2840 quad_iterator quad = triangulation.begin_quad(level),
2841 endc =
2842 (level == n_levels - 1 ?
2843 quad_iterator(triangulation.end_quad()) :
2844 triangulation.begin_quad(level + 1));
2845 for (; quad != endc; ++quad)
2846 {
2847 ++number_cache.n_quads_level[level];
2848 if (quad->has_children() == false)
2849 ++number_cache.n_active_quads_level[level];
2850 }
2851
2852 // update total number of quads
2853 number_cache.n_quads += number_cache.n_quads_level[level];
2854 number_cache.n_active_quads +=
2855 number_cache.n_active_quads_level[level];
2856 }
2857 }
2858 else
2859 {
2860 // for dim>2, there are no levels for quads
2861 number_cache.n_quads_level.clear();
2862 number_cache.n_active_quads_level.clear();
2863
2864 quad_iterator quad = triangulation.begin_quad(),
2865 endc = triangulation.end_quad();
2866 for (; quad != endc; ++quad)
2867 {
2868 ++number_cache.n_quads;
2869 if (quad->has_children() == false)
2870 ++number_cache.n_active_quads;
2871 }
2872 }
2873
2874 // wait for the background computation for lines
2875 update_lines.join();
2876 }
2877
2893 template <int dim, int spacedim>
2894 static void
2897 const unsigned int level_objects,
2899 {
2900 // update quads, lines and n_levels in number_cache. since we
2901 // don't access any of these numbers, we can do this in the
2902 // background
2903 Threads::Task<void> update_quads_and_lines = Threads::new_task(
2904 static_cast<
2905 void (*)(const Triangulation<dim, spacedim> &,
2906 const unsigned int,
2908 &compute_number_cache_dim<dim, spacedim>),
2910 level_objects,
2912 number_cache));
2913
2914 using hex_iterator =
2916
2917 //---------------------------------
2918 // update the number of hexes on the different levels in the
2919 // cache
2920 number_cache.n_hexes = 0;
2921 number_cache.n_active_hexes = 0;
2922
2923 // for 3d, hexes have levels so take count the objects per
2924 // level and globally
2925 if (dim == 3)
2926 {
2927 // count the number of levels; the function we called
2928 // above on a separate Task for quads (recursively, via
2929 // the lines function) also does this and puts it into
2930 // number_cache.n_levels, but this datum may not yet be
2931 // available as we call the function on a separate task
2932 unsigned int n_levels = 0;
2933 if (level_objects > 0)
2934 // find the last level on which there are used cells
2935 for (unsigned int level = 0; level < level_objects; ++level)
2936 if (triangulation.begin(level) != triangulation.end(level))
2937 n_levels = level + 1;
2938
2939 number_cache.n_hexes_level.resize(n_levels);
2940 number_cache.n_active_hexes_level.resize(n_levels);
2941
2942 for (unsigned int level = 0; level < n_levels; ++level)
2943 {
2944 // count hexes on this level
2945 number_cache.n_hexes_level[level] = 0;
2946 number_cache.n_active_hexes_level[level] = 0;
2947
2948 hex_iterator hex = triangulation.begin_hex(level),
2949 endc = (level == n_levels - 1 ?
2950 hex_iterator(triangulation.end_hex()) :
2951 triangulation.begin_hex(level + 1));
2952 for (; hex != endc; ++hex)
2953 {
2954 ++number_cache.n_hexes_level[level];
2955 if (hex->has_children() == false)
2956 ++number_cache.n_active_hexes_level[level];
2957 }
2958
2959 // update total number of hexes
2960 number_cache.n_hexes += number_cache.n_hexes_level[level];
2961 number_cache.n_active_hexes +=
2962 number_cache.n_active_hexes_level[level];
2963 }
2964 }
2965 else
2966 {
2967 // for dim>3, there are no levels for hexes
2968 number_cache.n_hexes_level.clear();
2969 number_cache.n_active_hexes_level.clear();
2970
2971 hex_iterator hex = triangulation.begin_hex(),
2972 endc = triangulation.end_hex();
2973 for (; hex != endc; ++hex)
2974 {
2975 ++number_cache.n_hexes;
2976 if (hex->has_children() == false)
2977 ++number_cache.n_active_hexes;
2978 }
2979 }
2980
2981 // wait for the background computation for quads
2982 update_quads_and_lines.join();
2983 }
2984
2985
2986 template <int dim, int spacedim>
2987 static void
2990 const unsigned int level_objects,
2992 {
2993 compute_number_cache_dim(triangulation, level_objects, number_cache);
2994
2995 number_cache.active_cell_index_partitioner =
2996 std::make_shared<const Utilities::MPI::Partitioner>(
2997 triangulation.n_active_cells());
2998
2999 number_cache.level_cell_index_partitioners.resize(
3000 triangulation.n_levels());
3001 for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
3002 number_cache.level_cell_index_partitioners[level] =
3003 std::make_shared<const Utilities::MPI::Partitioner>(
3004 triangulation.n_cells(level));
3005 }
3006
3007
3008 template <int spacedim>
3009 static void
3012
3013
3014 template <int dim, int spacedim>
3015 static void
3017 {
3018 // each face can be neighbored on two sides
3019 // by cells. according to the face's
3020 // intrinsic normal we define the left
3021 // neighbor as the one for which the face
3022 // normal points outward, and store that
3023 // one first; the second one is then
3024 // the right neighbor for which the
3025 // face normal points inward. This
3026 // information depends on the type of cell
3027 // and local number of face for the
3028 // 'standard ordering and orientation' of
3029 // faces and then on the face_orientation
3030 // information for the real mesh. Set up a
3031 // table to have fast access to those
3032 // offsets (0 for left and 1 for
3033 // right). Some of the values are invalid
3034 // as they reference too large face
3035 // numbers, but we just leave them at a
3036 // zero value.
3037 //
3038 // Note, that in 2d for lines as faces the
3039 // normal direction given in the
3040 // GeometryInfo class is not consistent. We
3041 // thus define here that the normal for a
3042 // line points to the right if the line
3043 // points upwards.
3044 //
3045 // There is one more point to
3046 // consider, however: if we have
3047 // dim<spacedim, then we may have
3048 // cases where cells are
3049 // inverted. In effect, both
3050 // cells think they are the left
3051 // neighbor of an edge, for
3052 // example, which leads us to
3053 // forget neighborship
3054 // information (a case that shows
3055 // this is
3056 // codim_one/hanging_nodes_02). We
3057 // store whether a cell is
3058 // inverted using the
3059 // direction_flag, so if a cell
3060 // has a false direction_flag,
3061 // then we need to invert our
3062 // selection whether we are a
3063 // left or right neighbor in all
3064 // following computations.
3065 //
3066 // first index: dimension (minus 2)
3067 // second index: local face index
3068 // third index: face_orientation (false and true)
3069 static const unsigned int left_right_offset[2][6][2] = {
3070 // quadrilateral
3071 {{0, 1}, // face 0, face_orientation = false and true
3072 {1, 0}, // face 1, face_orientation = false and true
3073 {1, 0}, // face 2, face_orientation = false and true
3074 {0, 1}, // face 3, face_orientation = false and true
3075 {0, 0}, // face 4, invalid face
3076 {0, 0}}, // face 5, invalid face
3077 // hexahedron
3078 {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
3079
3080 // now create a vector of the two active
3081 // neighbors (left and right) for each face
3082 // and fill it by looping over all cells. For
3083 // cases with anisotropic refinement and more
3084 // then one cell neighboring at a given side
3085 // of the face we will automatically get the
3086 // active one on the highest level as we loop
3087 // over cells from lower levels first.
3089 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
3090 adjacent_cells(2 * triangulation.n_raw_faces(), dummy);
3091
3092 for (const auto &cell : triangulation.cell_iterators())
3093 for (auto f : cell->face_indices())
3094 {
3096 cell->face(f);
3097
3098 const unsigned int offset =
3099 (cell->direction_flag() ?
3100 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
3101 1 -
3102 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
3103
3104 adjacent_cells[2 * face->index() + offset] = cell;
3105
3106 // if this cell is not refined, but the
3107 // face is, then we'll have to set our
3108 // cell as neighbor for the child faces
3109 // as well. Fortunately the normal
3110 // orientation of children will be just
3111 // the same.
3112 if (dim == 2)
3113 {
3114 if (cell->is_active() && face->has_children())
3115 {
3116 adjacent_cells[2 * face->child(0)->index() + offset] =
3117 cell;
3118 adjacent_cells[2 * face->child(1)->index() + offset] =
3119 cell;
3120 }
3121 }
3122 else // -> dim == 3
3123 {
3124 // We need the same as in 2d
3125 // here. Furthermore, if the face is
3126 // refined with cut_x or cut_y then
3127 // those children again in the other
3128 // direction, and if this cell is
3129 // refined isotropically (along the
3130 // face) then the neighbor will
3131 // (probably) be refined as cut_x or
3132 // cut_y along the face. For those
3133 // neighboring children cells, their
3134 // neighbor will be the current,
3135 // inactive cell, as our children are
3136 // too fine to be neighbors. Catch that
3137 // case by also acting on inactive
3138 // cells with isotropic refinement
3139 // along the face. If the situation
3140 // described is not present, the data
3141 // will be overwritten later on when we
3142 // visit cells on finer levels, so no
3143 // harm will be done.
3144 if (face->has_children() &&
3145 (cell->is_active() ||
3147 cell->refinement_case(), f) ==
3149 {
3150 for (unsigned int c = 0; c < face->n_children(); ++c)
3151 adjacent_cells[2 * face->child(c)->index() + offset] =
3152 cell;
3153 if (face->child(0)->has_children())
3154 {
3155 adjacent_cells[2 * face->child(0)->child(0)->index() +
3156 offset] = cell;
3157 adjacent_cells[2 * face->child(0)->child(1)->index() +
3158 offset] = cell;
3159 }
3160 if (face->child(1)->has_children())
3161 {
3162 adjacent_cells[2 * face->child(1)->child(0)->index() +
3163 offset] = cell;
3164 adjacent_cells[2 * face->child(1)->child(1)->index() +
3165 offset] = cell;
3166 }
3167 } // if cell active and face refined
3168 } // else -> dim==3
3169 } // for all faces of all cells
3170
3171 // now loop again over all cells and set the
3172 // corresponding neighbor cell. Note, that we
3173 // have to use the opposite of the
3174 // left_right_offset in this case as we want
3175 // the offset of the neighbor, not our own.
3176 for (const auto &cell : triangulation.cell_iterators())
3177 for (auto f : cell->face_indices())
3178 {
3179 const unsigned int offset =
3180 (cell->direction_flag() ?
3181 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
3182 1 -
3183 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
3184 cell->set_neighbor(
3185 f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]);
3186 }
3187 }
3188
3189
3193 template <int dim, int spacedim>
3194 static void
3196 const std::vector<CellData<dim>> &cells,
3197 const SubCellData &subcelldata,
3199 {
3200 AssertThrow(vertices.size() > 0, ExcMessage("No vertices given"));
3201 AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
3202
3203 // Check that all cells have positive volume.
3204#ifndef _MSC_VER
3205 // TODO: The following code does not compile with MSVC. Find a way
3206 // around it
3207 if (dim == spacedim)
3208 for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
3209 {
3210 // If we should check for distorted cells, then we permit them
3211 // to exist. If a cell has negative measure, then it must be
3212 // distorted (the converse is not necessarily true); hence
3213 // throw an exception if no such cells should exist.
3215 {
3216 const double cell_measure = GridTools::cell_measure<spacedim>(
3217 vertices,
3218 ArrayView<const unsigned int>(cells[cell_no].vertices));
3219 AssertThrow(cell_measure > 0, ExcGridHasInvalidCell(cell_no));
3220 }
3221 }
3222#endif
3223
3224 // clear old content
3225 tria.levels.clear();
3226 tria.levels.push_back(
3227 std::make_unique<
3229
3230 if (dim > 1)
3231 tria.faces = std::make_unique<
3233
3234 // copy vertices
3235 tria.vertices = vertices;
3236 tria.vertices_used.assign(vertices.size(), true);
3237
3238 // compute connectivity
3239 const auto connectivity = build_connectivity<unsigned int>(cells);
3240 const unsigned int n_cell = cells.size();
3241
3242 // TriaObjects: lines
3243 if (dim >= 2)
3244 {
3245 auto &lines_0 = tria.faces->lines; // data structure to be filled
3246
3247 // get connectivity between quads and lines
3248 const auto &crs = connectivity.entity_to_entities(1, 0);
3249 const unsigned int n_lines = crs.ptr.size() - 1;
3250
3251 // allocate memory
3252 reserve_space_(lines_0, n_lines);
3253
3254 // loop over lines
3255 for (unsigned int line = 0; line < n_lines; ++line)
3256 for (unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
3257 ++i, ++j)
3258 lines_0.cells[line * GeometryInfo<1>::faces_per_cell + j] =
3259 crs.col[i]; // set vertex indices
3260 }
3261
3262 // TriaObjects: quads
3263 if (dim == 3)
3264 {
3265 auto &quads_0 = tria.faces->quads; // data structures to be filled
3266 auto &faces = *tria.faces;
3267
3268 // get connectivity between quads and lines
3269 const auto &crs = connectivity.entity_to_entities(2, 1);
3270 const unsigned int n_quads = crs.ptr.size() - 1;
3271
3272 // allocate memory
3273 reserve_space_(quads_0, n_quads);
3274 reserve_space_(faces, 2 /*structdim*/, n_quads);
3275
3276 // loop over all quads -> entity type, line indices/orientations
3277 for (unsigned int q = 0, k = 0; q < n_quads; ++q)
3278 {
3279 // set entity type of quads
3280 faces.set_quad_type(q, connectivity.entity_types(2)[q]);
3281
3282 // loop over all its lines
3283 for (unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
3284 ++i, ++j, ++k)
3285 {
3286 // set line index
3287 quads_0.cells[q * GeometryInfo<3>::lines_per_face + j] =
3288 crs.col[i];
3289
3290 // set line orientations
3291 const unsigned char combined_orientation =
3292 connectivity.entity_orientations(1)
3293 .get_combined_orientation(k);
3294 // it doesn't make sense to set any flags except
3295 // orientation for a line
3296 Assert(
3297 combined_orientation ==
3299 combined_orientation ==
3302 faces.quads_line_orientations
3304 combined_orientation ==
3306 }
3307 }
3308 }
3309
3310 // TriaObjects/TriaLevel: cell
3311 {
3312 auto &cells_0 = tria.levels[0]->cells; // data structure to be filled
3313 auto &level = *tria.levels[0];
3314
3315 // get connectivity between cells/faces and cells/cells
3316 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
3317 const auto &nei = connectivity.entity_to_entities(dim, dim);
3318
3319 // in 2d optional: since in in pure QUAD meshes same line
3320 // orientations can be guaranteed
3321 bool orientation_needed = false;
3322 if (dim == 3)
3323 orientation_needed = true;
3324 else if (dim == 2)
3325 {
3326 const auto &orientations = connectivity.entity_orientations(1);
3327 for (unsigned int i = 0; i < orientations.n_objects(); ++i)
3328 if (orientations.get_combined_orientation(i) !=
3330 {
3331 orientation_needed = true;
3332 break;
3333 }
3334 }
3335
3336 // allocate memory
3337 reserve_space_(cells_0, n_cell);
3338 reserve_space_(level, spacedim, n_cell, orientation_needed);
3339
3340 // loop over all cells
3341 for (unsigned int cell = 0; cell < n_cell; ++cell)
3342 {
3343 // set material ids
3344 cells_0.boundary_or_material_id[cell].material_id =
3345 cells[cell].material_id;
3346
3347 // set manifold ids
3348 cells_0.manifold_id[cell] = cells[cell].manifold_id;
3349
3350 // set entity types
3351 level.reference_cell[cell] = connectivity.entity_types(dim)[cell];
3352
3353 // loop over faces
3354 for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
3355 ++i, ++j)
3356 {
3357 // set neighbor if not at boundary
3358 if (nei.col[i] != static_cast<unsigned int>(-1))
3359 level.neighbors[cell * GeometryInfo<dim>::faces_per_cell +
3360 j] = {0, nei.col[i]};
3361
3362 // set face indices
3363 cells_0.cells[cell * GeometryInfo<dim>::faces_per_cell + j] =
3364 crs.col[i];
3365
3366 // set face orientation if needed
3367 if (orientation_needed)
3368 {
3369 level.face_orientations.set_combined_orientation(
3371 connectivity.entity_orientations(dim - 1)
3372 .get_combined_orientation(i));
3373 }
3374 }
3375 }
3376 }
3377
3378 // TriaFaces: boundary id of boundary faces
3379 if (dim > 1)
3380 {
3381 auto &bids_face = dim == 3 ?
3382 tria.faces->quads.boundary_or_material_id :
3383 tria.faces->lines.boundary_or_material_id;
3384
3385 // count number of cells a face is belonging to
3386 std::vector<unsigned int> count(bids_face.size(), 0);
3387
3388 // get connectivity between cells/faces
3389 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
3390
3391 // count how many cells are adjacent to the same face
3392 for (unsigned int cell = 0; cell < cells.size(); ++cell)
3393 for (unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
3394 count[crs.col[i]]++;
3395
3396 // loop over all faces
3397 for (unsigned int face = 0; face < count.size(); ++face)
3398 {
3399 if (count[face] != 1) // inner face
3400 continue;
3401
3402 // boundary faces ...
3403 bids_face[face].boundary_id = 0;
3404
3405 if (dim != 3)
3406 continue;
3407
3408 // ... and the lines of quads in 3d
3409 const auto &crs = connectivity.entity_to_entities(2, 1);
3410 for (unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
3411 tria.faces->lines.boundary_or_material_id[crs.col[i]]
3412 .boundary_id = 0;
3413 }
3414 }
3415 else // 1d
3416 {
3417 static const unsigned int t_tba = static_cast<unsigned int>(-1);
3418 static const unsigned int t_inner = static_cast<unsigned int>(-2);
3419
3420 std::vector<unsigned int> type(vertices.size(), t_tba);
3421
3422 const auto &crs = connectivity.entity_to_entities(1, 0);
3423
3424 for (unsigned int cell = 0; cell < cells.size(); ++cell)
3425 for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
3426 ++i, ++j)
3427 if (type[crs.col[i]] != t_inner)
3428 type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
3429
3430 for (unsigned int face = 0; face < type.size(); ++face)
3431 {
3432 // note: we also treat manifolds here!?
3433 (*tria.vertex_to_manifold_id_map_1d)[face] =
3435 if (type[face] != t_inner && type[face] != t_tba)
3436 (*tria.vertex_to_boundary_id_map_1d)[face] = type[face];
3437 }
3438 }
3439
3440 // SubCellData: line
3441 if (dim >= 2)
3442 process_subcelldata(connectivity.entity_to_entities(1, 0),
3443 tria.faces->lines,
3444 subcelldata.boundary_lines,
3445 vertices);
3446
3447 // SubCellData: quad
3448 if (dim == 3)
3449 process_subcelldata(connectivity.entity_to_entities(2, 0),
3450 tria.faces->quads,
3451 subcelldata.boundary_quads,
3452 vertices);
3453 }
3454
3455
3456 template <int structdim, int spacedim, typename T>
3457 static void
3459 const CRS<T> &crs,
3460 TriaObjects &obj,
3461 const std::vector<CellData<structdim>> &boundary_objects_in,
3462 const std::vector<Point<spacedim>> &vertex_locations)
3463 {
3464 AssertDimension(obj.structdim, structdim);
3465
3466 if (boundary_objects_in.empty())
3467 return; // empty subcelldata -> nothing to do
3468
3469 // pre-sort subcelldata
3470 auto boundary_objects = boundary_objects_in;
3471
3472 // ... sort vertices
3473 for (auto &boundary_object : boundary_objects)
3474 std::sort(boundary_object.vertices.begin(),
3475 boundary_object.vertices.end());
3476
3477 // ... sort cells
3478 std::sort(boundary_objects.begin(),
3479 boundary_objects.end(),
3480 [](const auto &a, const auto &b) {
3481 return a.vertices < b.vertices;
3482 });
3483
3484 unsigned int counter = 0;
3485
3486 std::vector<unsigned int> key;
3488
3489 for (unsigned int o = 0; o < obj.n_objects(); ++o)
3490 {
3491 auto &boundary_id = obj.boundary_or_material_id[o].boundary_id;
3492 auto &manifold_id = obj.manifold_id[o];
3493
3494 // assert that object has not been visited yet and its value
3495 // has not been modified yet
3496 AssertThrow(boundary_id == 0 ||
3501
3502 // create key
3503 key.assign(crs.col.data() + crs.ptr[o],
3504 crs.col.data() + crs.ptr[o + 1]);
3505 std::sort(key.begin(), key.end());
3506
3507 // is subcelldata provided? -> binary search
3508 const auto subcell_object =
3509 std::lower_bound(boundary_objects.begin(),
3510 boundary_objects.end(),
3511 key,
3512 [&](const auto &cell, const auto &key) {
3513 return cell.vertices < key;
3514 });
3515
3516 // no subcelldata provided for this object
3517 if (subcell_object == boundary_objects.end() ||
3518 subcell_object->vertices != key)
3519 continue;
3520
3521 ++counter;
3522
3523 // set manifold id
3524 manifold_id = subcell_object->manifold_id;
3525
3526 // set boundary id
3527 if (subcell_object->boundary_id !=
3529 {
3530 (void)vertex_locations;
3533 ExcMessage(
3534 "The input arguments for creating a triangulation "
3535 "specified a boundary id for an internal face. This "
3536 "is not allowed."
3537 "\n\n"
3538 "The object in question has vertex indices " +
3539 [subcell_object]() {
3540 std::string s;
3541 for (const auto v : subcell_object->vertices)
3542 s += std::to_string(v) + ',';
3543 return s;
3544 }() +
3545 " which are located at positions " +
3546 [vertex_locations, subcell_object]() {
3547 std::ostringstream s;
3548 for (const auto v : subcell_object->vertices)
3549 s << '(' << vertex_locations[v] << ')';
3550 return s.str();
3551 }() +
3552 "."));
3553 boundary_id = subcell_object->boundary_id;
3554 }
3555 }
3556
3557 // make sure that all subcelldata entries have been processed
3558 // TODO: this is not guaranteed, why?
3559 // AssertDimension(counter, boundary_objects_in.size());
3560 (void)counter;
3561 }
3562
3563
3564
3565 static void
3567 const unsigned structdim,
3568 const unsigned int size)
3569 {
3570 const unsigned int dim = faces.dim;
3571
3572 const unsigned int max_lines_per_face = 2 * structdim;
3573
3574 if (dim == 3 && structdim == 2)
3575 {
3576 // quad entity types
3577 faces.quad_is_quadrilateral.assign(size, true);
3578
3579 // quad line orientations
3580 faces.quads_line_orientations.assign(size * max_lines_per_face,
3581 true);
3582 }
3583 }
3584
3585
3586
3587 static void
3589 const unsigned int spacedim,
3590 const unsigned int size,
3591 const bool orientation_needed)
3592 {
3593 const unsigned int dim = level.dim;
3594
3595 const unsigned int max_faces_per_cell = 2 * dim;
3596
3597 level.active_cell_indices.assign(size, numbers::invalid_unsigned_int);
3598 level.subdomain_ids.assign(size, 0);
3599 level.level_subdomain_ids.assign(size, 0);
3600
3601 level.refine_flags.assign(size, 0u);
3602 level.refine_choice.assign(size, 0u);
3603 level.coarsen_flags.assign(size, false);
3604
3605 level.parents.assign((size + 1) / 2, -1);
3606
3607 if (dim == spacedim - 1)
3608 level.direction_flags.assign(size, true);
3609
3610 level.neighbors.assign(size * max_faces_per_cell, {-1, -1});
3611
3612 level.reference_cell.assign(size, ReferenceCells::Invalid);
3613
3614 if (orientation_needed)
3615 level.face_orientations.reinit(size * max_faces_per_cell);
3616
3617
3618 level.global_active_cell_indices.assign(size,
3620 level.global_level_cell_indices.assign(size,
3622 }
3623
3624
3625
3626 static void
3627 reserve_space_(TriaObjects &obj, const unsigned int size)
3628 {
3629 const unsigned int structdim = obj.structdim;
3630
3631 const unsigned int max_children_per_cell = 1 << structdim;
3632 const unsigned int max_faces_per_cell = 2 * structdim;
3633
3634 obj.used.assign(size, true);
3635 obj.boundary_or_material_id.assign(
3636 size,
3638 BoundaryOrMaterialId());
3639 obj.manifold_id.assign(size, -1);
3640 obj.user_flags.assign(size, false);
3641 obj.user_data.resize(size);
3642
3643 if (structdim > 1) // TODO: why?
3644 obj.refinement_cases.assign(size, 0);
3645
3646 obj.children.assign(max_children_per_cell / 2 * size, -1);
3647
3648 obj.cells.assign(max_faces_per_cell * size, -1);
3649
3650 if (structdim <= 2)
3651 {
3652 obj.next_free_single = size - 1;
3653 obj.next_free_pair = 0;
3655 }
3656 else
3657 {
3658 obj.next_free_single = obj.next_free_pair = 0;
3659 }
3660 }
3661
3662
3678 template <int spacedim>
3679 static void
3682 std::vector<unsigned int> &,
3683 std::vector<unsigned int> &)
3684 {
3685 const unsigned int dim = 1;
3686
3687 // first we need to reset the
3688 // neighbor pointers of the
3689 // neighbors of this cell's
3690 // children to this cell. This is
3691 // different for one dimension,
3692 // since there neighbors can have a
3693 // refinement level differing from
3694 // that of this cell's children by
3695 // more than one level.
3696
3697 Assert(!cell->child(0)->has_children() &&
3698 !cell->child(1)->has_children(),
3700
3701 // first do it for the cells to the
3702 // left
3703 if (cell->neighbor(0).state() == IteratorState::valid)
3704 if (cell->neighbor(0)->has_children())
3705 {
3707 cell->neighbor(0);
3708 Assert(neighbor->level() == cell->level(), ExcInternalError());
3709
3710 // right child
3711 neighbor = neighbor->child(1);
3712 while (true)
3713 {
3714 Assert(neighbor->neighbor(1) == cell->child(0),
3716 neighbor->set_neighbor(1, cell);
3717
3718 // move on to further
3719 // children on the
3720 // boundary between this
3721 // cell and its neighbor
3722 if (neighbor->has_children())
3723 neighbor = neighbor->child(1);
3724 else
3725 break;
3726 }
3727 }
3728
3729 // now do it for the cells to the
3730 // left
3731 if (cell->neighbor(1).state() == IteratorState::valid)
3732 if (cell->neighbor(1)->has_children())
3733 {
3735 cell->neighbor(1);
3736 Assert(neighbor->level() == cell->level(), ExcInternalError());
3737
3738 // left child
3739 neighbor = neighbor->child(0);
3740 while (true)
3741 {
3742 Assert(neighbor->neighbor(0) == cell->child(1),
3744 neighbor->set_neighbor(0, cell);
3745
3746 // move on to further
3747 // children on the
3748 // boundary between this
3749 // cell and its neighbor
3750 if (neighbor->has_children())
3751 neighbor = neighbor->child(0);
3752 else
3753 break;
3754 }
3755 }
3756
3757
3758 // delete the vertex which will not
3759 // be needed anymore. This vertex
3760 // is the second of the first child
3761 triangulation.vertices_used[cell->child(0)->vertex_index(1)] = false;
3762
3763 // invalidate children. clear user
3764 // pointers, to avoid that they may
3765 // appear at unwanted places later
3766 // on...
3767 for (unsigned int child = 0; child < cell->n_children(); ++child)
3768 {
3769 cell->child(child)->clear_user_data();
3770 cell->child(child)->clear_user_flag();
3771 cell->child(child)->clear_used_flag();
3772 }
3773
3774
3775 // delete pointer to children
3776 cell->clear_children();
3777 cell->clear_user_flag();
3778 }
3779
3780
3781
3782 template <int spacedim>
3783 static void
3786 std::vector<unsigned int> &line_cell_count,
3787 std::vector<unsigned int> &)
3788 {
3789 const unsigned int dim = 2;
3790 const RefinementCase<dim> ref_case = cell->refinement_case();
3791
3792 Assert(line_cell_count.size() == triangulation.n_raw_lines(),
3794
3795 // vectors to hold all lines which
3796 // may be deleted
3797 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3798 lines_to_delete(0);
3799
3800 lines_to_delete.reserve(4 * 2 + 4);
3801
3802 // now we decrease the counters for
3803 // lines contained in the child
3804 // cells
3805 for (unsigned int c = 0; c < cell->n_children(); ++c)
3806 {
3808 cell->child(c);
3809 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3810 --line_cell_count[child->line_index(l)];
3811 }
3812
3813
3814 // delete the vertex which will not
3815 // be needed anymore. This vertex
3816 // is the second of the second line
3817 // of the first child, if the cell
3818 // is refined with cut_xy, else there
3819 // is no inner vertex.
3820 // additionally delete unneeded inner
3821 // lines
3822 if (ref_case == RefinementCase<dim>::cut_xy)
3823 {
3825 .vertices_used[cell->child(0)->line(1)->vertex_index(1)] = false;
3826
3827 lines_to_delete.push_back(cell->child(0)->line(1));
3828 lines_to_delete.push_back(cell->child(0)->line(3));
3829 lines_to_delete.push_back(cell->child(3)->line(0));
3830 lines_to_delete.push_back(cell->child(3)->line(2));
3831 }
3832 else
3833 {
3834 unsigned int inner_face_no =
3835 ref_case == RefinementCase<dim>::cut_x ? 1 : 3;
3836
3837 // the inner line will not be
3838 // used any more
3839 lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
3840 }
3841
3842 // invalidate children
3843 for (unsigned int child = 0; child < cell->n_children(); ++child)
3844 {
3845 cell->child(child)->clear_user_data();
3846 cell->child(child)->clear_user_flag();
3847 cell->child(child)->clear_used_flag();
3848 }
3849
3850
3851 // delete pointer to children
3852 cell->clear_children();
3853 cell->clear_refinement_case();
3854 cell->clear_user_flag();
3855
3856 // look at the refinement of outer
3857 // lines. if nobody needs those
3858 // anymore we can add them to the
3859 // list of lines to be deleted.
3860 for (unsigned int line_no = 0;
3861 line_no < GeometryInfo<dim>::lines_per_cell;
3862 ++line_no)
3863 {
3865 cell->line(line_no);
3866
3867 if (line->has_children())
3868 {
3869 // if one of the cell counters is
3870 // zero, the other has to be as well
3871
3872 Assert((line_cell_count[line->child_index(0)] == 0 &&
3873 line_cell_count[line->child_index(1)] == 0) ||
3874 (line_cell_count[line->child_index(0)] > 0 &&
3875 line_cell_count[line->child_index(1)] > 0),
3877
3878 if (line_cell_count[line->child_index(0)] == 0)
3879 {
3880 for (unsigned int c = 0; c < 2; ++c)
3881 Assert(!line->child(c)->has_children(),
3883
3884 // we may delete the line's
3885 // children and the middle vertex
3886 // as no cell references them
3887 // anymore
3889 .vertices_used[line->child(0)->vertex_index(1)] = false;
3890
3891 lines_to_delete.push_back(line->child(0));
3892 lines_to_delete.push_back(line->child(1));
3893
3894 line->clear_children();
3895 }
3896 }
3897 }
3898
3899 // finally, delete unneeded lines
3900
3901 // clear user pointers, to avoid that
3902 // they may appear at unwanted places
3903 // later on...
3904 // same for user flags, then finally
3905 // delete the lines
3906 typename std::vector<
3908 line = lines_to_delete.begin(),
3909 endline = lines_to_delete.end();
3910 for (; line != endline; ++line)
3911 {
3912 (*line)->clear_user_data();
3913 (*line)->clear_user_flag();
3914 (*line)->clear_used_flag();
3915 }
3916 }
3917
3918
3919
3920 template <int spacedim>
3921 static void
3924 std::vector<unsigned int> &line_cell_count,
3925 std::vector<unsigned int> &quad_cell_count)
3926 {
3927 const unsigned int dim = 3;
3928
3929 Assert(line_cell_count.size() == triangulation.n_raw_lines(),
3931 Assert(quad_cell_count.size() == triangulation.n_raw_quads(),
3933
3934 // first of all, we store the RefineCase of
3935 // this cell
3936 const RefinementCase<dim> ref_case = cell->refinement_case();
3937 // vectors to hold all lines and quads which
3938 // may be deleted
3939 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3940 lines_to_delete(0);
3941 std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3942 quads_to_delete(0);
3943
3944 lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3945 quads_to_delete.reserve(6 * 4 + 12);
3946
3947 // now we decrease the counters for lines and
3948 // quads contained in the child cells
3949 for (unsigned int c = 0; c < cell->n_children(); ++c)
3950 {
3952 cell->child(c);
3953 const auto line_indices = TriaAccessorImplementation::
3954 Implementation::get_line_indices_of_cell(*child);
3955 for (const unsigned int l : cell->line_indices())
3956 --line_cell_count[line_indices[l]];
3957 for (auto f : GeometryInfo<dim>::face_indices())
3958 --quad_cell_count[child->quad_index(f)];
3959 }
3960
3961 //-------------------------------------
3962 // delete interior quads and lines and the
3963 // interior vertex, depending on the
3964 // refinement case of the cell
3965 //
3966 // for append quads and lines: only append
3967 // them to the list of objects to be deleted
3968
3969 switch (ref_case)
3970 {
3972 quads_to_delete.push_back(cell->child(0)->face(1));
3973 break;
3975 quads_to_delete.push_back(cell->child(0)->face(3));
3976 break;
3978 quads_to_delete.push_back(cell->child(0)->face(5));
3979 break;
3981 quads_to_delete.push_back(cell->child(0)->face(1));
3982 quads_to_delete.push_back(cell->child(0)->face(3));
3983 quads_to_delete.push_back(cell->child(3)->face(0));
3984 quads_to_delete.push_back(cell->child(3)->face(2));
3985
3986 lines_to_delete.push_back(cell->child(0)->line(11));
3987 break;
3989 quads_to_delete.push_back(cell->child(0)->face(1));
3990 quads_to_delete.push_back(cell->child(0)->face(5));
3991 quads_to_delete.push_back(cell->child(3)->face(0));
3992 quads_to_delete.push_back(cell->child(3)->face(4));
3993
3994 lines_to_delete.push_back(cell->child(0)->line(5));
3995 break;
3997 quads_to_delete.push_back(cell->child(0)->face(3));
3998 quads_to_delete.push_back(cell->child(0)->face(5));
3999 quads_to_delete.push_back(cell->child(3)->face(2));
4000 quads_to_delete.push_back(cell->child(3)->face(4));
4001
4002 lines_to_delete.push_back(cell->child(0)->line(7));
4003 break;
4005 quads_to_delete.push_back(cell->child(0)->face(1));
4006 quads_to_delete.push_back(cell->child(2)->face(1));
4007 quads_to_delete.push_back(cell->child(4)->face(1));
4008 quads_to_delete.push_back(cell->child(6)->face(1));
4009
4010 quads_to_delete.push_back(cell->child(0)->face(3));
4011 quads_to_delete.push_back(cell->child(1)->face(3));
4012 quads_to_delete.push_back(cell->child(4)->face(3));
4013 quads_to_delete.push_back(cell->child(5)->face(3));
4014
4015 quads_to_delete.push_back(cell->child(0)->face(5));
4016 quads_to_delete.push_back(cell->child(1)->face(5));
4017 quads_to_delete.push_back(cell->child(2)->face(5));
4018 quads_to_delete.push_back(cell->child(3)->face(5));
4019
4020 lines_to_delete.push_back(cell->child(0)->line(5));
4021 lines_to_delete.push_back(cell->child(0)->line(7));
4022 lines_to_delete.push_back(cell->child(0)->line(11));
4023 lines_to_delete.push_back(cell->child(7)->line(0));
4024 lines_to_delete.push_back(cell->child(7)->line(2));
4025 lines_to_delete.push_back(cell->child(7)->line(8));
4026 // delete the vertex which will not
4027 // be needed anymore. This vertex
4028 // is the vertex at the heart of
4029 // this cell, which is the sixth of
4030 // the first child
4031 triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
4032 false;
4033 break;
4034 default:
4035 // only remaining case is
4036 // no_refinement, thus an error
4038 break;
4039 }
4040
4041
4042 // invalidate children
4043 for (unsigned int child = 0; child < cell->n_children(); ++child)
4044 {
4045 cell->child(child)->clear_user_data();
4046 cell->child(child)->clear_user_flag();
4047
4048 for (auto f : GeometryInfo<dim>::face_indices())
4049 // set flags denoting deviations from standard orientation of
4050 // faces back to initialization values
4051 cell->child(child)->set_combined_face_orientation(
4053
4054 cell->child(child)->clear_used_flag();
4055 }
4056
4057
4058 // delete pointer to children
4059 cell->clear_children();
4060 cell->clear_refinement_case();
4061 cell->clear_user_flag();
4062
4063 // so far we only looked at inner quads,
4064 // lines and vertices. Now we have to
4065 // consider outer ones as well. here, we have
4066 // to check, whether there are other cells
4067 // still needing these objects. otherwise we
4068 // can delete them. first for quads (and
4069 // their inner lines).
4070
4071 for (const unsigned int quad_no : GeometryInfo<dim>::face_indices())
4072 {
4074 cell->face(quad_no);
4075
4076 Assert(
4077 (GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) &&
4078 quad->has_children()) ||
4079 GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) ==
4082
4083 switch (quad->refinement_case())
4084 {
4085 case RefinementCase<dim - 1>::no_refinement:
4086 // nothing to do as the quad
4087 // is not refined
4088 break;
4089 case RefinementCase<dim - 1>::cut_x:
4090 case RefinementCase<dim - 1>::cut_y:
4091 {
4092 // if one of the cell counters is
4093 // zero, the other has to be as
4094 // well
4095 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
4096 quad_cell_count[quad->child_index(1)] == 0) ||
4097 (quad_cell_count[quad->child_index(0)] > 0 &&
4098 quad_cell_count[quad->child_index(1)] > 0),
4100 // it might be, that the quad is
4101 // refined twice anisotropically,
4102 // first check, whether we may
4103 // delete possible grand_children
4104 unsigned int deleted_grandchildren = 0;
4105 unsigned int number_of_child_refinements = 0;
4106
4107 for (unsigned int c = 0; c < 2; ++c)
4108 if (quad->child(c)->has_children())
4109 {
4110 ++number_of_child_refinements;
4111 // if one of the cell counters is
4112 // zero, the other has to be as
4113 // well
4114 Assert(
4115 (quad_cell_count[quad->child(c)->child_index(0)] ==
4116 0 &&
4117 quad_cell_count[quad->child(c)->child_index(1)] ==
4118 0) ||
4119 (quad_cell_count[quad->child(c)->child_index(0)] >
4120 0 &&
4121 quad_cell_count[quad->child(c)->child_index(1)] >
4122 0),
4124 if (quad_cell_count[quad->child(c)->child_index(0)] ==
4125 0)
4126 {
4127 // Assert, that the two
4128 // anisotropic
4129 // refinements add up to
4130 // isotropic refinement
4131 Assert(quad->refinement_case() +
4132 quad->child(c)->refinement_case() ==
4135 // we may delete the
4136 // quad's children and
4137 // the inner line as no
4138 // cell references them
4139 // anymore
4140 quads_to_delete.push_back(
4141 quad->child(c)->child(0));
4142 quads_to_delete.push_back(
4143 quad->child(c)->child(1));
4144 if (quad->child(c)->refinement_case() ==
4146 lines_to_delete.push_back(
4147 quad->child(c)->child(0)->line(1));
4148 else
4149 lines_to_delete.push_back(
4150 quad->child(c)->child(0)->line(3));
4151 quad->child(c)->clear_children();
4152 quad->child(c)->clear_refinement_case();
4153 ++deleted_grandchildren;
4154 }
4155 }
4156 // if no grandchildren are left, we
4157 // may as well delete the
4158 // refinement of the inner line
4159 // between our children and the
4160 // corresponding vertex
4161 if (number_of_child_refinements > 0 &&
4162 deleted_grandchildren == number_of_child_refinements)
4163 {
4165 middle_line;
4166 if (quad->refinement_case() == RefinementCase<2>::cut_x)
4167 middle_line = quad->child(0)->line(1);
4168 else
4169 middle_line = quad->child(0)->line(3);
4170
4171 lines_to_delete.push_back(middle_line->child(0));
4172 lines_to_delete.push_back(middle_line->child(1));
4174 .vertices_used[middle_vertex_index<dim, spacedim>(
4175 middle_line)] = false;
4176 middle_line->clear_children();
4177 }
4178
4179 // now consider the direct children
4180 // of the given quad
4181 if (quad_cell_count[quad->child_index(0)] == 0)
4182 {
4183 // we may delete the quad's
4184 // children and the inner line
4185 // as no cell references them
4186 // anymore
4187 quads_to_delete.push_back(quad->child(0));
4188 quads_to_delete.push_back(quad->child(1));
4189 if (quad->refinement_case() == RefinementCase<2>::cut_x)
4190 lines_to_delete.push_back(quad->child(0)->line(1));
4191 else
4192 lines_to_delete.push_back(quad->child(0)->line(3));
4193
4194 // if the counters just dropped
4195 // to zero, otherwise the
4196 // children would have been
4197 // deleted earlier, then this
4198 // cell's children must have
4199 // contained the anisotropic
4200 // quad children. thus, if
4201 // those have again anisotropic
4202 // children, which are in
4203 // effect isotropic children of
4204 // the original quad, those are
4205 // still needed by a
4206 // neighboring cell and we
4207 // cannot delete them. instead,
4208 // we have to reset this quad's
4209 // refine case to isotropic and
4210 // set the children
4211 // accordingly.
4212 if (quad->child(0)->has_children())
4213 if (quad->refinement_case() ==
4215 {
4216 // now evereything is
4217 // quite complicated. we
4218 // have the children
4219 // numbered according to
4220 //
4221 // *---*---*
4222 // |n+1|m+1|
4223 // *---*---*
4224 // | n | m |
4225 // *---*---*
4226 //
4227 // from the original
4228 // anisotropic
4229 // refinement. we have to
4230 // reorder them as
4231 //
4232 // *---*---*
4233 // | m |m+1|
4234 // *---*---*
4235 // | n |n+1|
4236 // *---*---*
4237 //
4238 // for isotropic refinement.
4239 //
4240 // this is a bit ugly, of
4241 // course: loop over all
4242 // cells on all levels
4243 // and look for faces n+1
4244 // (switch_1) and m
4245 // (switch_2).
4246 const typename Triangulation<dim, spacedim>::
4247 quad_iterator switch_1 =
4248 quad->child(0)->child(1),
4249 switch_2 =
4250 quad->child(1)->child(0);
4251
4252 Assert(!switch_1->has_children(),
4254 Assert(!switch_2->has_children(),
4256
4257 const int switch_1_index = switch_1->index();
4258 const int switch_2_index = switch_2->index();
4259 for (unsigned int l = 0;
4260 l < triangulation.levels.size();
4261 ++l)
4262 for (unsigned int h = 0;
4263 h <
4264 triangulation.levels[l]->cells.n_objects();
4265 ++h)
4266 for (const unsigned int q :
4268 {
4269 const int index =
4271 ->cells.get_bounding_object_indices(
4272 h)[q];
4273 if (index == switch_1_index)
4274 triangulation.levels[l]
4275 ->cells.get_bounding_object_indices(
4276 h)[q] = switch_2_index;
4277 else if (index == switch_2_index)
4278 triangulation.levels[l]
4279 ->cells.get_bounding_object_indices(
4280 h)[q] = switch_1_index;
4281 }
4282 // now we have to copy
4283 // all information of the
4284 // two quads
4285 const int switch_1_lines[4] = {
4286 static_cast<signed int>(
4287 switch_1->line_index(0)),
4288 static_cast<signed int>(
4289 switch_1->line_index(1)),
4290 static_cast<signed int>(
4291 switch_1->line_index(2)),
4292 static_cast<signed int>(
4293 switch_1->line_index(3))};
4294 const bool switch_1_line_orientations[4] = {
4295 switch_1->line_orientation(0),
4296 switch_1->line_orientation(1),
4297 switch_1->line_orientation(2),
4298 switch_1->line_orientation(3)};
4299 const types::boundary_id switch_1_boundary_id =
4300 switch_1->boundary_id();
4301 const unsigned int switch_1_user_index =
4302 switch_1->user_index();
4303 const bool switch_1_user_flag =
4304 switch_1->user_flag_set();
4305
4306 switch_1->set_bounding_object_indices(
4307 {switch_2->line_index(0),
4308 switch_2->line_index(1),
4309 switch_2->line_index(2),
4310 switch_2->line_index(3)});
4311 switch_1->set_line_orientation(
4312 0, switch_2->line_orientation(0));
4313 switch_1->set_line_orientation(
4314 1, switch_2->line_orientation(1));
4315 switch_1->set_line_orientation(
4316 2, switch_2->line_orientation(2));
4317 switch_1->set_line_orientation(
4318 3, switch_2->line_orientation(3));
4319 switch_1->set_boundary_id_internal(
4320 switch_2->boundary_id());
4321 switch_1->set_manifold_id(
4322 switch_2->manifold_id());
4323 switch_1->set_user_index(switch_2->user_index());
4324 if (switch_2->user_flag_set())
4325 switch_1->set_user_flag();
4326 else
4327 switch_1->clear_user_flag();
4328
4329 switch_2->set_bounding_object_indices(
4330 {switch_1_lines[0],
4331 switch_1_lines[1],
4332 switch_1_lines[2],
4333 switch_1_lines[3]});
4334 switch_2->set_line_orientation(
4335 0, switch_1_line_orientations[0]);
4336 switch_2->set_line_orientation(
4337 1, switch_1_line_orientations[1]);
4338 switch_2->set_line_orientation(
4339 2, switch_1_line_orientations[2]);
4340 switch_2->set_line_orientation(
4341 3, switch_1_line_orientations[3]);
4342 switch_2->set_boundary_id_internal(
4343 switch_1_boundary_id);
4344 switch_2->set_manifold_id(
4345 switch_1->manifold_id());
4346 switch_2->set_user_index(switch_1_user_index);
4347 if (switch_1_user_flag)
4348 switch_2->set_user_flag();
4349 else
4350 switch_2->clear_user_flag();
4351
4352 const unsigned int child_0 =
4353 quad->child(0)->child_index(0);
4354 const unsigned int child_2 =
4355 quad->child(1)->child_index(0);
4356 quad->clear_children();
4357 quad->clear_refinement_case();
4358 quad->set_refinement_case(
4360 quad->set_children(0, child_0);
4361 quad->set_children(2, child_2);
4362 std::swap(quad_cell_count[child_0 + 1],
4363 quad_cell_count[child_2]);
4364 }
4365 else
4366 {
4367 // the face was refined
4368 // with cut_y, thus the
4369 // children are already
4370 // in correct order. we
4371 // only have to set them
4372 // correctly, deleting
4373 // the indirection of two
4374 // anisotropic refinement
4375 // and going directly
4376 // from the quad to
4377 // isotropic children
4378 const unsigned int child_0 =
4379 quad->child(0)->child_index(0);
4380 const unsigned int child_2 =
4381 quad->child(1)->child_index(0);
4382 quad->clear_children();
4383 quad->clear_refinement_case();
4384 quad->set_refinement_case(
4386 quad->set_children(0, child_0);
4387 quad->set_children(2, child_2);
4388 }
4389 else
4390 {
4391 quad->clear_children();
4392 quad->clear_refinement_case();
4393 }
4394 }
4395 break;
4396 }
4397 case RefinementCase<dim - 1>::cut_xy:
4398 {
4399 // if one of the cell counters is
4400 // zero, the others have to be as
4401 // well
4402
4403 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
4404 quad_cell_count[quad->child_index(1)] == 0 &&
4405 quad_cell_count[quad->child_index(2)] == 0 &&
4406 quad_cell_count[quad->child_index(3)] == 0) ||
4407 (quad_cell_count[quad->child_index(0)] > 0 &&
4408 quad_cell_count[quad->child_index(1)] > 0 &&
4409 quad_cell_count[quad->child_index(2)] > 0 &&
4410 quad_cell_count[quad->child_index(3)] > 0),
4412
4413 if (quad_cell_count[quad->child_index(0)] == 0)
4414 {
4415 // we may delete the quad's
4416 // children, the inner lines
4417 // and the middle vertex as no
4418 // cell references them anymore
4419 lines_to_delete.push_back(quad->child(0)->line(1));
4420 lines_to_delete.push_back(quad->child(3)->line(0));
4421 lines_to_delete.push_back(quad->child(0)->line(3));
4422 lines_to_delete.push_back(quad->child(3)->line(2));
4423
4424 for (unsigned int child = 0; child < quad->n_children();
4425 ++child)
4426 quads_to_delete.push_back(quad->child(child));
4427
4429 .vertices_used[quad->child(0)->vertex_index(3)] =
4430 false;
4431
4432 quad->clear_children();
4433 quad->clear_refinement_case();
4434 }
4435 }
4436 break;
4437
4438 default:
4440 break;
4441 }
4442 }
4443
4444 // now we repeat a similar procedure
4445 // for the outer lines of this cell.
4446
4447 // if in debug mode: check that each
4448 // of the lines for which we consider
4449 // deleting the children in fact has
4450 // children (the bits/coarsening_3d
4451 // test tripped over this initially)
4452 for (unsigned int line_no = 0;
4453 line_no < GeometryInfo<dim>::lines_per_cell;
4454 ++line_no)
4455 {
4457 cell->line(line_no);
4458
4459 Assert(
4460 (GeometryInfo<dim>::line_refinement_case(ref_case, line_no) &&
4461 line->has_children()) ||
4462 GeometryInfo<dim>::line_refinement_case(ref_case, line_no) ==
4465
4466 if (line->has_children())
4467 {
4468 // if one of the cell counters is
4469 // zero, the other has to be as well
4470
4471 Assert((line_cell_count[line->child_index(0)] == 0 &&
4472 line_cell_count[line->child_index(1)] == 0) ||
4473 (line_cell_count[line->child_index(0)] > 0 &&
4474 line_cell_count[line->child_index(1)] > 0),
4476
4477 if (line_cell_count[line->child_index(0)] == 0)
4478 {
4479 for (unsigned int c = 0; c < 2; ++c)
4480 Assert(!line->child(c)->has_children(),
4482
4483 // we may delete the line's
4484 // children and the middle vertex
4485 // as no cell references them
4486 // anymore
4488 .vertices_used[line->child(0)->vertex_index(1)] = false;
4489
4490 lines_to_delete.push_back(line->child(0));
4491 lines_to_delete.push_back(line->child(1));
4492
4493 line->clear_children();
4494 }
4495 }
4496 }
4497
4498 // finally, delete unneeded quads and lines
4499
4500 // clear user pointers, to avoid that
4501 // they may appear at unwanted places
4502 // later on...
4503 // same for user flags, then finally
4504 // delete the quads and lines
4505 typename std::vector<
4507 line = lines_to_delete.begin(),
4508 endline = lines_to_delete.end();
4509 for (; line != endline; ++line)
4510 {
4511 (*line)->clear_user_data();
4512 (*line)->clear_user_flag();
4513 (*line)->clear_used_flag();
4514 }
4515
4516 typename std::vector<
4518 quad = quads_to_delete.begin(),
4519 endquad = quads_to_delete.end();
4520 for (; quad != endquad; ++quad)
4521 {
4522 (*quad)->clear_user_data();
4523 (*quad)->clear_children();
4524 (*quad)->clear_refinement_case();
4525 (*quad)->clear_user_flag();
4526 (*quad)->clear_used_flag();
4527 }
4528 }
4529
4530
4548 template <int spacedim>
4549 static void
4552 unsigned int &next_unused_vertex,
4554 &next_unused_line,
4556 &next_unused_cell,
4557 const typename Triangulation<2, spacedim>::cell_iterator &cell)
4558 {
4559 const unsigned int dim = 2;
4560 // clear refinement flag
4561 const RefinementCase<dim> ref_case = cell->refine_flag_set();
4562 cell->clear_refine_flag();
4563
4564 /* For the refinement process: since we go the levels up from the
4565 lowest, there are (unlike above) only two possibilities: a neighbor
4566 cell is on the same level or one level up (in both cases, it may or
4567 may not be refined later on, but we don't care here).
4568
4569 First:
4570 Set up an array of the 3x3 vertices, which are distributed on the
4571 cell (the array consists of indices into the @p{vertices} std::vector
4572
4573 2--7--3
4574 | | |
4575 4--8--5
4576 | | |
4577 0--6--1
4578
4579 note: in case of cut_x or cut_y not all these vertices are needed for
4580 the new cells
4581
4582 Second:
4583 Set up an array of the new lines (the array consists of iterator
4584 pointers into the lines arrays)
4585
4586 .-6-.-7-. The directions are: .->-.->-.
4587 1 9 3 ^ ^ ^
4588 .-10.11-. .->-.->-.
4589 0 8 2 ^ ^ ^
4590 .-4-.-5-. .->-.->-.
4591
4592 cut_x:
4593 .-4-.-5-.
4594 | | |
4595 0 6 1
4596 | | |
4597 .-2-.-3-.
4598
4599 cut_y:
4600 .---5---.
4601 1 3
4602 .---6---.
4603 0 2
4604 .---4---.
4605
4606
4607 Third:
4608 Set up an array of neighbors:
4609
4610 6 7
4611 .--.--.
4612 1| | |3
4613 .--.--.
4614 0| | |2
4615 .--.--.
4616 4 5
4617
4618 We need this array for two reasons: first to get the lines which will
4619 bound the four subcells (if the neighboring cell is refined, these
4620 lines already exist), and second to update neighborship information.
4621 Since if a neighbor is not refined, its neighborship record only
4622 points to the present, unrefined, cell rather than the children we
4623 are presently creating, we only need the neighborship information
4624 if the neighbor cells are refined. In all other cases, we store
4625 the unrefined neighbor address
4626
4627 We also need for every neighbor (if refined) which number among its
4628 neighbors the present (unrefined) cell has, since that number is to
4629 be replaced and because that also is the number of the subline which
4630 will be the interface between that neighbor and the to be created
4631 cell. We will store this number (between 0 and 3) in the field
4632 @p{neighbors_neighbor}.
4633
4634 It would be sufficient to use the children of the common line to the
4635 neighbor, if we only wanted to get the new sublines and the new
4636 vertex, but because we need to update the neighborship information of
4637 the two refined subcells of the neighbor, we need to search these
4638 anyway.
4639
4640 Convention:
4641 The created children are numbered like this:
4642
4643 .--.--.
4644 |2 . 3|
4645 .--.--.
4646 |0 | 1|
4647 .--.--.
4648 */
4649 // collect the indices of the eight surrounding vertices
4650 // 2--7--3
4651 // | | |
4652 // 4--8--5
4653 // | | |
4654 // 0--6--1
4655 int new_vertices[9];
4656 for (unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
4657 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
4658 for (unsigned int line_no = 0; line_no < 4; ++line_no)
4659 if (cell->line(line_no)->has_children())
4660 new_vertices[4 + line_no] =
4661 cell->line(line_no)->child(0)->vertex_index(1);
4662
4663 if (ref_case == RefinementCase<dim>::cut_xy)
4664 {
4665 // find the next
4666 // unused vertex and
4667 // allocate it for
4668 // the new vertex we
4669 // need here
4670 while (triangulation.vertices_used[next_unused_vertex] == true)
4671 ++next_unused_vertex;
4672 Assert(next_unused_vertex < triangulation.vertices.size(),
4673 ExcMessage(
4674 "Internal error: During refinement, the triangulation "
4675 "wants to access an element of the 'vertices' array "
4676 "but it turns out that the array is not large enough."));
4677 triangulation.vertices_used[next_unused_vertex] = true;
4678
4679 new_vertices[8] = next_unused_vertex;
4680
4681 // determine middle vertex by transfinite interpolation to be
4682 // consistent with what happens to quads in a
4683 // Triangulation<3,3> when they are refined
4684 triangulation.vertices[next_unused_vertex] =
4685 cell->center(true, true);
4686 }
4687
4688
4689 // Now the lines:
4691 unsigned int lmin = 8;
4692 unsigned int lmax = 12;
4693 if (ref_case != RefinementCase<dim>::cut_xy)
4694 {
4695 lmin = 6;
4696 lmax = 7;
4697 }
4698
4699 for (unsigned int l = lmin; l < lmax; ++l)
4700 {
4701 while (next_unused_line->used() == true)
4702 ++next_unused_line;
4703 new_lines[l] = next_unused_line;
4704 ++next_unused_line;
4705
4706 AssertIsNotUsed(new_lines[l]);
4707 }
4708
4709 if (ref_case == RefinementCase<dim>::cut_xy)
4710 {
4711 // .-6-.-7-.
4712 // 1 9 3
4713 // .-10.11-.
4714 // 0 8 2
4715 // .-4-.-5-.
4716
4717 // lines 0-7 already exist, create only the four interior
4718 // lines 8-11
4719 unsigned int l = 0;
4720 for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
4721 for (unsigned int c = 0; c < 2; ++c, ++l)
4722 new_lines[l] = cell->line(face_no)->child(c);
4723 Assert(l == 8, ExcInternalError());
4724
4725 new_lines[8]->set_bounding_object_indices(
4726 {new_vertices[6], new_vertices[8]});
4727 new_lines[9]->set_bounding_object_indices(
4728 {new_vertices[8], new_vertices[7]});
4729 new_lines[10]->set_bounding_object_indices(
4730 {new_vertices[4], new_vertices[8]});
4731 new_lines[11]->set_bounding_object_indices(
4732 {new_vertices[8], new_vertices[5]});
4733 }
4734 else if (ref_case == RefinementCase<dim>::cut_x)
4735 {
4736 // .-4-.-5-.
4737 // | | |
4738 // 0 6 1
4739 // | | |
4740 // .-2-.-3-.
4741 new_lines[0] = cell->line(0);
4742 new_lines[1] = cell->line(1);
4743 new_lines[2] = cell->line(2)->child(0);
4744 new_lines[3] = cell->line(2)->child(1);
4745 new_lines[4] = cell->line(3)->child(0);
4746 new_lines[5] = cell->line(3)->child(1);
4747 new_lines[6]->set_bounding_object_indices(
4748 {new_vertices[6], new_vertices[7]});
4749 }
4750 else
4751 {
4753 // .---5---.
4754 // 1 3
4755 // .---6---.
4756 // 0 2
4757 // .---4---.
4758 new_lines[0] = cell->line(0)->child(0);
4759 new_lines[1] = cell->line(0)->child(1);
4760 new_lines[2] = cell->line(1)->child(0);
4761 new_lines[3] = cell->line(1)->child(1);
4762 new_lines[4] = cell->line(2);
4763 new_lines[5] = cell->line(3);
4764 new_lines[6]->set_bounding_object_indices(
4765 {new_vertices[4], new_vertices[5]});
4766 }
4767
4768 for (unsigned int l = lmin; l < lmax; ++l)
4769 {
4770 new_lines[l]->set_used_flag();
4771 new_lines[l]->clear_user_flag();
4772 new_lines[l]->clear_user_data();
4773 new_lines[l]->clear_children();
4774 // interior line
4775 new_lines[l]->set_boundary_id_internal(
4777 new_lines[l]->set_manifold_id(cell->manifold_id());
4778 }
4779
4780 // Now add the four (two)
4781 // new cells!
4784 while (next_unused_cell->used() == true)
4785 ++next_unused_cell;
4786
4787 const unsigned int n_children = GeometryInfo<dim>::n_children(ref_case);
4788 for (unsigned int i = 0; i < n_children; ++i)
4789 {
4790 AssertIsNotUsed(next_unused_cell);
4791 subcells[i] = next_unused_cell;
4792 ++next_unused_cell;
4793 if (i % 2 == 1 && i < n_children - 1)
4794 while (next_unused_cell->used() == true)
4795 ++next_unused_cell;
4796 }
4797
4798 if (ref_case == RefinementCase<dim>::cut_xy)
4799 {
4800 // children:
4801 // .--.--.
4802 // |2 . 3|
4803 // .--.--.
4804 // |0 | 1|
4805 // .--.--.
4806 // lines:
4807 // .-6-.-7-.
4808 // 1 9 3
4809 // .-10.11-.
4810 // 0 8 2
4811 // .-4-.-5-.
4812 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4813 new_lines[8]->index(),
4814 new_lines[4]->index(),
4815 new_lines[10]->index()});
4816 subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
4817 new_lines[2]->index(),
4818 new_lines[5]->index(),
4819 new_lines[11]->index()});
4820 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
4821 new_lines[9]->index(),
4822 new_lines[10]->index(),
4823 new_lines[6]->index()});
4824 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
4825 new_lines[3]->index(),
4826 new_lines[11]->index(),
4827 new_lines[7]->index()});
4828 }
4829 else if (ref_case == RefinementCase<dim>::cut_x)
4830 {
4831 // children:
4832 // .--.--.
4833 // | . |
4834 // .0 . 1.
4835 // | | |
4836 // .--.--.
4837 // lines:
4838 // .-4-.-5-.
4839 // | | |
4840 // 0 6 1
4841 // | | |
4842 // .-2-.-3-.
4843 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4844 new_lines[6]->index(),
4845 new_lines[2]->index(),
4846 new_lines[4]->index()});
4847 subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
4848 new_lines[1]->index(),
4849 new_lines[3]->index(),
4850 new_lines[5]->index()});
4851 }
4852 else
4853 {
4855 // children:
4856 // .-----.
4857 // | 1 |
4858 // .-----.
4859 // | 0 |
4860 // .-----.
4861 // lines:
4862 // .---5---.
4863 // 1 3
4864 // .---6---.
4865 // 0 2
4866 // .---4---.
4867 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4868 new_lines[2]->index(),
4869 new_lines[4]->index(),
4870 new_lines[6]->index()});
4871 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
4872 new_lines[3]->index(),
4873 new_lines[6]->index(),
4874 new_lines[5]->index()});
4875 }
4876
4877 types::subdomain_id subdomainid = cell->subdomain_id();
4878
4879 for (unsigned int i = 0; i < n_children; ++i)
4880 {
4881 subcells[i]->set_used_flag();
4882 subcells[i]->clear_refine_flag();
4883 subcells[i]->clear_user_flag();
4884 subcells[i]->clear_user_data();
4885 subcells[i]->clear_children();
4886 // inherit material properties
4887 subcells[i]->set_material_id(cell->material_id());
4888 subcells[i]->set_manifold_id(cell->manifold_id());
4889 subcells[i]->set_subdomain_id(subdomainid);
4890
4891 if (i % 2 == 0)
4892 subcells[i]->set_parent(cell->index());
4893 }
4894
4895
4896
4897 // set child index for even children i=0,2 (0)
4898 for (unsigned int i = 0; i < n_children / 2; ++i)
4899 cell->set_children(2 * i, subcells[2 * i]->index());
4900 // set the refine case
4901 cell->set_refinement_case(ref_case);
4902
4903 // note that the
4904 // refinement flag was
4905 // already cleared at the
4906 // beginning of this function
4907
4908 if (dim == spacedim - 1)
4909 for (unsigned int c = 0; c < n_children; ++c)
4910 cell->child(c)->set_direction_flag(cell->direction_flag());
4911 }
4912
4913
4914
4915 template <int dim, int spacedim>
4918 const bool check_for_distorted_cells)
4919 {
4920 AssertDimension(dim, 2);
4921
4922 // Check whether a new level is needed. We have to check for
4923 // this on the highest level only
4924 for (const auto &cell : triangulation.active_cell_iterators_on_level(
4925 triangulation.levels.size() - 1))
4926 if (cell->refine_flag_set())
4927 {
4928 triangulation.levels.push_back(
4929 std::make_unique<
4931 break;
4932 }
4933
4936 line != triangulation.end_line();
4937 ++line)
4938 {
4939 line->clear_user_flag();
4940 line->clear_user_data();
4941 }
4942
4943 unsigned int n_single_lines = 0;
4944 unsigned int n_lines_in_pairs = 0;
4945 unsigned int needed_vertices = 0;
4946
4947 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4948 {
4949 // count number of flagged cells on this level and compute
4950 // how many new vertices and new lines will be needed
4951 unsigned int needed_cells = 0;
4952
4953 for (const auto &cell :
4954 triangulation.active_cell_iterators_on_level(level))
4955 if (cell->refine_flag_set())
4956 {
4957 if (cell->reference_cell() == ReferenceCells::Triangle)
4958 {
4959 needed_cells += 4;
4960 needed_vertices += 0;
4961 n_single_lines += 3;
4962 }
4963 else if (cell->reference_cell() ==
4965 {
4966 needed_cells += 4;
4967 needed_vertices += 1;
4968 n_single_lines += 4;
4969 }
4970 else
4971 {
4973 }
4974
4975 for (const auto line_no : cell->face_indices())
4976 {
4977 auto line = cell->line(line_no);
4978 if (line->has_children() == false)
4979 line->set_user_flag();
4980 }
4981 }
4982
4983
4984 const unsigned int used_cells =
4985 std::count(triangulation.levels[level + 1]->cells.used.begin(),
4986 triangulation.levels[level + 1]->cells.used.end(),
4987 true);
4988
4989
4990 reserve_space(*triangulation.levels[level + 1],
4991 used_cells + needed_cells,
4992 2,
4993 spacedim);
4994
4995 reserve_space(triangulation.levels[level + 1]->cells,
4996 needed_cells,
4997 0);
4998 }
4999
5000 for (auto line = triangulation.begin_line();
5001 line != triangulation.end_line();
5002 ++line)
5003 if (line->user_flag_set())
5004 {
5005 Assert(line->has_children() == false, ExcInternalError());
5006 n_lines_in_pairs += 2;
5007 needed_vertices += 1;
5008 }
5009
5010 reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
5011
5012 needed_vertices += std::count(triangulation.vertices_used.begin(),
5013 triangulation.vertices_used.end(),
5014 true);
5015
5016 if (needed_vertices > triangulation.vertices.size())
5017 {
5018 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
5019 triangulation.vertices_used.resize(needed_vertices, false);
5020 }
5021
5022 unsigned int next_unused_vertex = 0;
5023
5024 {
5027 endl = triangulation.end_line();
5029 next_unused_line = triangulation.begin_raw_line();
5030
5031 for (; line != endl; ++line)
5032 if (line->user_flag_set())
5033 {
5034 // This line needs to be refined. Find the next unused vertex
5035 // and set it appropriately
5036 while (triangulation.vertices_used[next_unused_vertex] == true)
5037 ++next_unused_vertex;
5038 Assert(next_unused_vertex < triangulation.vertices.size(),
5039 ExcMessage(
5040 "Internal error: During refinement, the triangulation "
5041 "wants to access an element of the 'vertices' array "
5042 "but it turns out that the array is not large "
5043 "enough."));
5044 triangulation.vertices_used[next_unused_vertex] = true;
5045
5046 triangulation.vertices[next_unused_vertex] = line->center(true);
5047
5048 bool pair_found = false;
5049 (void)pair_found;
5050 for (; next_unused_line != endl; ++next_unused_line)
5051 if (!next_unused_line->used() &&
5052 !(++next_unused_line)->used())
5053 {
5054 --next_unused_line;
5055 pair_found = true;
5056 break;
5057 }
5058 Assert(pair_found, ExcInternalError());
5059
5060 line->set_children(0, next_unused_line->index());
5061
5063 children[2] = {next_unused_line, ++next_unused_line};
5064
5065 AssertIsNotUsed(children[0]);
5066 AssertIsNotUsed(children[1]);
5067
5068 children[0]->set_bounding_object_indices(
5069 {line->vertex_index(0), next_unused_vertex});
5070 children[1]->set_bounding_object_indices(
5071 {next_unused_vertex, line->vertex_index(1)});
5072
5073 for (auto &child : children)
5074 {
5075 child->set_used_flag();
5076 child->clear_children();
5077 child->clear_user_data();
5078 child->clear_user_flag();
5079 child->set_boundary_id_internal(line->boundary_id());
5080 child->set_manifold_id(line->manifold_id());
5081 // Line orientation is relative to the cell it is on so
5082 // those cannot be set at this point.
5083 }
5084
5085 line->clear_user_flag();
5086 }
5087 }
5088
5089 reserve_space(triangulation.faces->lines, 0, n_single_lines);
5090
5092 cells_with_distorted_children;
5093
5095 next_unused_line = triangulation.begin_raw_line();
5096
5097 const auto create_children = [](auto &triangulation,
5098 unsigned int &next_unused_vertex,
5099 auto &next_unused_line,
5100 auto &next_unused_cell,
5101 const auto &cell) {
5102 const auto ref_case = cell->refine_flag_set();
5103 cell->clear_refine_flag();
5104
5105 unsigned int n_new_vertices = 0;
5106
5107 if (cell->reference_cell() == ReferenceCells::Triangle)
5108 n_new_vertices = 6;
5109 else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
5110 n_new_vertices = 9;
5111 else
5113
5114 std::vector<unsigned int> new_vertices(n_new_vertices,
5116 for (unsigned int vertex_no = 0; vertex_no < cell->n_vertices();
5117 ++vertex_no)
5118 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
5119 for (unsigned int line_no = 0; line_no < cell->n_lines(); ++line_no)
5120 if (cell->line(line_no)->has_children())
5121 new_vertices[cell->n_vertices() + line_no] =
5122 cell->line(line_no)->child(0)->vertex_index(1);
5123
5124 if (cell->reference_cell() == ReferenceCells::Quadrilateral)
5125 {
5126 while (triangulation.vertices_used[next_unused_vertex] == true)
5127 ++next_unused_vertex;
5128 Assert(
5129 next_unused_vertex < triangulation.vertices.size(),
5130 ExcMessage(
5131 "Internal error: During refinement, the triangulation wants "
5132 "to access an element of the 'vertices' array but it turns "
5133 "out that the array is not large enough."));
5134 triangulation.vertices_used[next_unused_vertex] = true;
5135
5136 new_vertices[8] = next_unused_vertex;
5137
5138 triangulation.vertices[next_unused_vertex] =
5139 cell->center(true, true);
5140 }
5141
5142 std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
5143 12>
5144 new_lines;
5145 std::array<unsigned char, 12> inherited_orientations;
5146 inherited_orientations.fill(
5148 unsigned int lmin = 0;
5149 unsigned int lmax = 0;
5150
5151 if (cell->reference_cell() == ReferenceCells::Triangle)
5152 {
5153 lmin = 6;
5154 lmax = 9;
5155 // For triangles, the innermost faces are always reversed for the
5156 // first three children and are in the standard orientation for
5157 // the last one.
5158 std::fill(inherited_orientations.begin() + lmin,
5159 inherited_orientations.begin() + lmax,
5161 }
5162 else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
5163 {
5164 lmin = 8;
5165 lmax = 12;
5166 }
5167 else
5168 {
5170 }
5171
5172 for (unsigned int l = lmin; l < lmax; ++l)
5173 {
5174 while (next_unused_line->used() == true)
5175 ++next_unused_line;
5176 new_lines[l] = next_unused_line;
5177 ++next_unused_line;
5178
5179 AssertIsNotUsed(new_lines[l]);
5180 }
5181
5182 // set up lines which have parents:
5183 for (const unsigned int face_no : cell->face_indices())
5184 {
5185 // Check the face (line) orientation to ensure that the (six or
5186 // eight) outer lines in new_lines are indexed in the default
5187 // orientation. This way we can index into this array in the
5188 // without special casing orientations (e.g., quadrilateral child
5189 // 3 will always have lines 9, 3, 11, 7) when setting child lines.
5190 const unsigned char combined_orientation =
5191 cell->combined_face_orientation(face_no);
5192 Assert(combined_orientation ==
5194 combined_orientation ==
5197 for (unsigned int c = 0; c < 2; ++c)
5198 {
5199 new_lines[2 * face_no + c] = cell->line(face_no)->child(c);
5200 inherited_orientations[2 * face_no + c] =
5201 cell->combined_face_orientation(face_no);
5202 }
5203 if (combined_orientation ==
5205 std::swap(new_lines[2 * face_no], new_lines[2 * face_no + 1]);
5206 }
5207
5208 // set up lines which do not have parents:
5209 if (cell->reference_cell() == ReferenceCells::Triangle)
5210 {
5211 new_lines[6]->set_bounding_object_indices(
5212 {new_vertices[3], new_vertices[4]});
5213 new_lines[7]->set_bounding_object_indices(
5214 {new_vertices[4], new_vertices[5]});
5215 new_lines[8]->set_bounding_object_indices(
5216 {new_vertices[5], new_vertices[3]});
5217 }
5218 else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
5219 {
5220 new_lines[8]->set_bounding_object_indices(
5221 {new_vertices[6], new_vertices[8]});
5222 new_lines[9]->set_bounding_object_indices(
5223 {new_vertices[8], new_vertices[7]});
5224 new_lines[10]->set_bounding_object_indices(
5225 {new_vertices[4], new_vertices[8]});
5226 new_lines[11]->set_bounding_object_indices(
5227 {new_vertices[8], new_vertices[5]});
5228 }
5229 else
5230 {
5232 }
5233
5234 for (unsigned int l = lmin; l < lmax; ++l)
5235 {
5236 new_lines[l]->set_used_flag();
5237 new_lines[l]->clear_user_flag();
5238 new_lines[l]->clear_user_data();
5239 new_lines[l]->clear_children();
5240 // new lines are always internal.
5241 new_lines[l]->set_boundary_id_internal(
5243 new_lines[l]->set_manifold_id(cell->manifold_id());
5244 }
5245
5248 while (next_unused_cell->used() == true)
5249 ++next_unused_cell;
5250
5251 unsigned int n_children = 0;
5252 if (cell->reference_cell() == ReferenceCells::Triangle)
5253 n_children = 4;
5254 else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
5255 n_children = 4;
5256 else
5258
5259 for (unsigned int i = 0; i < n_children; ++i)
5260 {
5261 AssertIsNotUsed(next_unused_cell);
5262 subcells[i] = next_unused_cell;
5263 ++next_unused_cell;
5264 if (i % 2 == 1 && i < n_children - 1)
5265 while (next_unused_cell->used() == true)
5266 ++next_unused_cell;
5267 }
5268
5269 // Assign lines to child cells:
5270 constexpr unsigned int X = numbers::invalid_unsigned_int;
5271 static constexpr ::ndarray<unsigned int, 4, 4> tri_child_lines =
5272 {{{{0, 8, 5, X}}, {{1, 2, 6, X}}, {{7, 3, 4, X}}, {{6, 7, 8, X}}}};
5273 static constexpr ::ndarray<unsigned int, 4, 4>
5274 quad_child_lines = {{{{0, 8, 4, 10}},
5275 {{8, 2, 5, 11}},
5276 {{1, 9, 10, 6}},
5277 {{9, 3, 11, 7}}}};
5278 // Here and below we assume that child cells have the same reference
5279 // cell type as the parent.
5280 const auto &child_lines =
5281 cell->reference_cell() == ReferenceCells::Triangle ?
5282 tri_child_lines :
5283 quad_child_lines;
5284 for (unsigned int i = 0; i < n_children; ++i)
5285 {
5286 if (cell->reference_cell() == ReferenceCells::Triangle)
5287 subcells[i]->set_bounding_object_indices(
5288 {new_lines[child_lines[i][0]]->index(),
5289 new_lines[child_lines[i][1]]->index(),
5290 new_lines[child_lines[i][2]]->index()});
5291 else
5292 subcells[i]->set_bounding_object_indices(
5293 {new_lines[child_lines[i][0]]->index(),
5294 new_lines[child_lines[i][1]]->index(),
5295 new_lines[child_lines[i][2]]->index(),
5296 new_lines[child_lines[i][3]]->index()});
5297
5298 subcells[i]->set_used_flag();
5299 subcells[i]->clear_refine_flag();
5300 subcells[i]->clear_user_flag();
5301 subcells[i]->clear_user_data();
5302 subcells[i]->clear_children();
5303 // inherit material properties
5304 subcells[i]->set_material_id(cell->material_id());
5305 subcells[i]->set_manifold_id(cell->manifold_id());
5306 subcells[i]->set_subdomain_id(cell->subdomain_id());
5307
5308 triangulation.levels[subcells[i]->level()]
5309 ->reference_cell[subcells[i]->index()] = cell->reference_cell();
5310
5311 // Finally, now that children are marked as used, we can set
5312 // orientation flags:
5313 for (unsigned int face_no : cell->face_indices())
5314 subcells[i]->set_combined_face_orientation(
5315 face_no, inherited_orientations[child_lines[i][face_no]]);
5316
5317 if (i % 2 == 0)
5318 subcells[i]->set_parent(cell->index());
5319 }
5320
5321 // Unlike the same lines on other children, the innermost triangle's
5322 // faces are all in the default orientation:
5323 if (cell->reference_cell() == ReferenceCells::Triangle)
5324 for (unsigned int face_no : cell->face_indices())
5325 subcells[3]->set_combined_face_orientation(
5327
5328 for (unsigned int i = 0; i < n_children / 2; ++i)
5329 cell->set_children(2 * i, subcells[2 * i]->index());
5330
5331 cell->set_refinement_case(ref_case);
5332
5333 if (dim == spacedim - 1)
5334 for (unsigned int c = 0; c < n_children; ++c)
5335 cell->child(c)->set_direction_flag(cell->direction_flag());
5336 };
5337
5338 for (int level = 0;
5339 level < static_cast<int>(triangulation.levels.size()) - 1;
5340 ++level)
5341 {
5343 next_unused_cell = triangulation.begin_raw(level + 1);
5344
5345 for (const auto &cell :
5346 triangulation.active_cell_iterators_on_level(level))
5347 if (cell->refine_flag_set())
5348 {
5349 create_children(triangulation,
5350 next_unused_vertex,
5351 next_unused_line,
5352 next_unused_cell,
5353 cell);
5354
5355 if (cell->reference_cell() == ReferenceCells::Quadrilateral &&
5356 check_for_distorted_cells &&
5357 has_distorted_children<dim, spacedim>(cell))
5358 cells_with_distorted_children.distorted_cells.push_back(
5359 cell);
5360
5361 triangulation.signals.post_refinement_on_cell(cell);
5362 }
5363 }
5364
5365 return cells_with_distorted_children;
5366 }
5367
5368
5369
5374 template <int spacedim>
5377 const bool /*check_for_distorted_cells*/)
5378 {
5379 const unsigned int dim = 1;
5380
5381 // Check whether a new level is needed. We have to check for
5382 // this on the highest level only
5383 for (const auto &cell : triangulation.active_cell_iterators_on_level(
5384 triangulation.levels.size() - 1))
5385 if (cell->refine_flag_set())
5386 {
5387 triangulation.levels.push_back(
5388 std::make_unique<
5390 break;
5391 }
5392
5393
5394 // check how much space is needed on every level. We need not
5395 // check the highest level since either - on the highest level
5396 // no cells are flagged for refinement - there are, but
5397 // prepare_refinement added another empty level
5398 unsigned int needed_vertices = 0;
5399 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
5400 {
5401 // count number of flagged
5402 // cells on this level
5403 unsigned int flagged_cells = 0;
5404
5405 for (const auto &acell :
5406 triangulation.active_cell_iterators_on_level(level))
5407 if (acell->refine_flag_set())
5408 ++flagged_cells;
5409
5410 // count number of used cells
5411 // on the next higher level
5412 const unsigned int used_cells =
5413 std::count(triangulation.levels[level + 1]->cells.used.begin(),
5414 triangulation.levels[level + 1]->cells.used.end(),
5415 true);
5416
5417 // reserve space for the used_cells cells already existing
5418 // on the next higher level as well as for the
5419 // 2*flagged_cells that will be created on that level
5420 reserve_space(*triangulation.levels[level + 1],
5422 flagged_cells,
5423 1,
5424 spacedim);
5425 // reserve space for 2*flagged_cells new lines on the next
5426 // higher level
5427 reserve_space(triangulation.levels[level + 1]->cells,
5429 flagged_cells,
5430 0);
5431
5432 needed_vertices += flagged_cells;
5433 }
5434
5435 // add to needed vertices how many
5436 // vertices are already in use
5437 needed_vertices += std::count(triangulation.vertices_used.begin(),
5438 triangulation.vertices_used.end(),
5439 true);
5440 // if we need more vertices: create them, if not: leave the
5441 // array as is, since shrinking is not really possible because
5442 // some of the vertices at the end may be in use
5443 if (needed_vertices > triangulation.vertices.size())
5444 {
5445 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
5446 triangulation.vertices_used.resize(needed_vertices, false);
5447 }
5448
5449
5450 // Do REFINEMENT on every level; exclude highest level as
5451 // above
5452
5453 // index of next unused vertex
5454 unsigned int next_unused_vertex = 0;
5455
5456 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
5457 {
5459 next_unused_cell = triangulation.begin_raw(level + 1);
5460
5461 for (const auto &cell :
5462 triangulation.active_cell_iterators_on_level(level))
5463 if (cell->refine_flag_set())
5464 {
5465 // clear refinement flag
5466 cell->clear_refine_flag();
5467
5468 // search for next unused
5469 // vertex
5470 while (triangulation.vertices_used[next_unused_vertex] ==
5471 true)
5472 ++next_unused_vertex;
5473 Assert(
5474 next_unused_vertex < triangulation.vertices.size(),
5475 ExcMessage(
5476 "Internal error: During refinement, the triangulation "
5477 "wants to access an element of the 'vertices' array "
5478 "but it turns out that the array is not large enough."));
5479
5480 // Now we always ask the cell itself where to put
5481 // the new point. The cell in turn will query the
5482 // manifold object internally.
5483 triangulation.vertices[next_unused_vertex] =
5484 cell->center(true);
5485
5486 triangulation.vertices_used[next_unused_vertex] = true;
5487
5488 // search for next two unused cell (++ takes care of
5489 // the end of the vector)
5491 first_child,
5492 second_child;
5493 while (next_unused_cell->used() == true)
5494 ++next_unused_cell;
5495 first_child = next_unused_cell;
5496 first_child->set_used_flag();
5497 first_child->clear_user_data();
5498 ++next_unused_cell;
5499 AssertIsNotUsed(next_unused_cell);
5500 second_child = next_unused_cell;
5501 second_child->set_used_flag();
5502 second_child->clear_user_data();
5503
5504 types::subdomain_id subdomainid = cell->subdomain_id();
5505
5506 // insert first child
5507 cell->set_children(0, first_child->index());
5508 first_child->clear_children();
5509 first_child->set_bounding_object_indices(
5510 {cell->vertex_index(0), next_unused_vertex});
5511 first_child->set_material_id(cell->material_id());
5512 first_child->set_manifold_id(cell->manifold_id());
5513 first_child->set_subdomain_id(subdomainid);
5514 if (dim == spacedim - 1)
5515 first_child->set_direction_flag(cell->direction_flag());
5516
5517 first_child->set_parent(cell->index());
5518
5519 // Set manifold id of the right face. Only do this
5520 // on the first child.
5521 first_child->face(1)->set_manifold_id(cell->manifold_id());
5522
5523 // reset neighborship info (refer to
5524 // internal::TriangulationImplementation::TriaLevel<0> for
5525 // details)
5526 first_child->set_neighbor(1, second_child);
5527 if (cell->neighbor(0).state() != IteratorState::valid)
5528 first_child->set_neighbor(0, cell->neighbor(0));
5529 else if (cell->neighbor(0)->is_active())
5530 {
5531 // since the neighbors level is always <=level,
5532 // if the cell is active, then there are no
5533 // cells to the left which may want to know
5534 // about this new child cell.
5535 Assert(cell->neighbor(0)->level() <= cell->level(),
5537 first_child->set_neighbor(0, cell->neighbor(0));
5538 }
5539 else
5540 // left neighbor is refined
5541 {
5542 // set neighbor to cell on same level
5543 const unsigned int nbnb = cell->neighbor_of_neighbor(0);
5544 first_child->set_neighbor(0,
5545 cell->neighbor(0)->child(nbnb));
5546
5547 // reset neighbor info of all right descendant
5548 // of the left neighbor of cell
5550 left_neighbor = cell->neighbor(0);
5551 while (left_neighbor->has_children())
5552 {
5553 left_neighbor = left_neighbor->child(nbnb);
5554 left_neighbor->set_neighbor(nbnb, first_child);
5555 }
5556 }
5557
5558 // insert second child
5559 second_child->clear_children();
5560 second_child->set_bounding_object_indices(
5561 {next_unused_vertex, cell->vertex_index(1)});
5562 second_child->set_neighbor(0, first_child);
5563 second_child->set_material_id(cell->material_id());
5564 second_child->set_manifold_id(cell->manifold_id());
5565 second_child->set_subdomain_id(subdomainid);
5566 if (dim == spacedim - 1)
5567 second_child->set_direction_flag(cell->direction_flag());
5568
5569 if (cell->neighbor(1).state() != IteratorState::valid)
5570 second_child->set_neighbor(1, cell->neighbor(1));
5571 else if (cell->neighbor(1)->is_active())
5572 {
5573 Assert(cell->neighbor(1)->level() <= cell->level(),
5575 second_child->set_neighbor(1, cell->neighbor(1));
5576 }
5577 else
5578 // right neighbor is refined same as above
5579 {
5580 const unsigned int nbnb = cell->neighbor_of_neighbor(1);
5581 second_child->set_neighbor(
5582 1, cell->neighbor(1)->child(nbnb));
5583
5585 right_neighbor = cell->neighbor(1);
5586 while (right_neighbor->has_children())
5587 {
5588 right_neighbor = right_neighbor->child(nbnb);
5589 right_neighbor->set_neighbor(nbnb, second_child);
5590 }
5591 }
5592 // inform all listeners that cell refinement is done
5593 triangulation.signals.post_refinement_on_cell(cell);
5594 }
5595 }
5596
5597 // in 1d, we can not have distorted children unless the parent
5598 // was already distorted (that is because we don't use
5599 // boundary information for 1d triangulations). so return an
5600 // empty list
5602 }
5603
5604
5609 template <int spacedim>
5612 const bool check_for_distorted_cells)
5613 {
5614 const unsigned int dim = 2;
5615
5616 // First check whether we can get away with isotropic refinement, or
5617 // whether we need to run through the full anisotropic algorithm
5618 {
5619 bool do_isotropic_refinement = true;
5620 for (const auto &cell : triangulation.active_cell_iterators())
5621 if (cell->refine_flag_set() == RefinementCase<dim>::cut_x ||
5622 cell->refine_flag_set() == RefinementCase<dim>::cut_y)
5623 {
5624 do_isotropic_refinement = false;
5625 break;
5626 }
5627
5628 if (do_isotropic_refinement)
5629 return execute_refinement_isotropic(triangulation,
5630 check_for_distorted_cells);
5631 }
5632
5633 // If we get here, we are doing anisotropic refinement.
5634
5635 // Check whether a new level is needed. We have to check for
5636 // this on the highest level only
5637 for (const auto &cell : triangulation.active_cell_iterators_on_level(
5638 triangulation.levels.size() - 1))
5639 if (cell->refine_flag_set())
5640 {
5641 triangulation.levels.push_back(
5642 std::make_unique<
5644 break;
5645 }
5646
5647 // TODO[WB]: we clear user flags and pointers of lines; we're going
5648 // to use them to flag which lines need refinement
5651 line != triangulation.end_line();
5652 ++line)
5653 {
5654 line->clear_user_flag();
5655 line->clear_user_data();
5656 }
5657 // running over all cells and lines count the number
5658 // n_single_lines of lines which can be stored as single
5659 // lines, e.g. inner lines
5660 unsigned int n_single_lines = 0;
5661
5662 // New lines to be created: number lines which are stored in
5663 // pairs (the children of lines must be stored in pairs)
5664 unsigned int n_lines_in_pairs = 0;
5665
5666 // check how much space is needed on every level. We need not
5667 // check the highest level since either - on the highest level
5668 // no cells are flagged for refinement - there are, but
5669 // prepare_refinement added another empty level
5670 unsigned int needed_vertices = 0;
5671 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
5672 {
5673 // count number of flagged cells on this level and compute
5674 // how many new vertices and new lines will be needed
5675 unsigned int needed_cells = 0;
5676
5677 for (const auto &cell :
5678 triangulation.active_cell_iterators_on_level(level))
5679 if (cell->refine_flag_set())
5680 {
5681 if (cell->refine_flag_set() == RefinementCase<dim>::cut_xy)
5682 {
5683 needed_cells += 4;
5684
5685 // new vertex at center of cell is needed in any
5686 // case
5687 ++needed_vertices;
5688
5689 // the four inner lines can be stored as singles
5690 n_single_lines += 4;
5691 }
5692 else // cut_x || cut_y
5693 {
5694 // set the flag showing that anisotropic
5695 // refinement is used for at least one cell
5696 triangulation.anisotropic_refinement = true;
5697
5698 needed_cells += 2;
5699 // no vertex at center
5700
5701 // the inner line can be stored as single
5702 n_single_lines += 1;
5703 }
5704
5705 // mark all faces (lines) for refinement; checking
5706 // locally whether the neighbor would also like to
5707 // refine them is rather difficult for lines so we
5708 // only flag them and after visiting all cells, we
5709 // decide which lines need refinement;
5710 for (const unsigned int line_no :
5712 {
5714 cell->refine_flag_set(), line_no) ==
5716 {
5718 line = cell->line(line_no);
5719 if (line->has_children() == false)
5720 line->set_user_flag();
5721 }
5722 }
5723 }
5724
5725
5726 // count number of used cells on the next higher level
5727 const unsigned int used_cells =
5728 std::count(triangulation.levels[level + 1]->cells.used.begin(),
5729 triangulation.levels[level + 1]->cells.used.end(),
5730 true);
5731
5732
5733 // reserve space for the used_cells cells already existing
5734 // on the next higher level as well as for the
5735 // needed_cells that will be created on that level
5736 reserve_space(*triangulation.levels[level + 1],
5737 used_cells + needed_cells,
5738 2,
5739 spacedim);
5740
5741 // reserve space for needed_cells new quads on the next
5742 // higher level
5743 reserve_space(triangulation.levels[level + 1]->cells,
5744 needed_cells,
5745 0);
5746 }
5747
5748 // now count the lines which were flagged for refinement
5751 line != triangulation.end_line();
5752 ++line)
5753 if (line->user_flag_set())
5754 {
5755 Assert(line->has_children() == false, ExcInternalError());
5756 n_lines_in_pairs += 2;
5757 needed_vertices += 1;
5758 }
5759 // reserve space for n_lines_in_pairs new lines. note, that
5760 // we can't reserve space for the single lines here as well,
5761 // as all the space reserved for lines in pairs would be
5762 // counted as unused and we would end up with too little space
5763 // to store all lines. memory reservation for n_single_lines
5764 // can only be done AFTER we refined the lines of the current
5765 // cells
5766 reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
5767
5768 // add to needed vertices how many vertices are already in use
5769 needed_vertices += std::count(triangulation.vertices_used.begin(),
5770 triangulation.vertices_used.end(),
5771 true);
5772 // if we need more vertices: create them, if not: leave the
5773 // array as is, since shrinking is not really possible because
5774 // some of the vertices at the end may be in use
5775 if (needed_vertices > triangulation.vertices.size())
5776 {
5777 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
5778 triangulation.vertices_used.resize(needed_vertices, false);
5779 }
5780
5781
5782 // Do REFINEMENT on every level; exclude highest level as
5783 // above
5784
5785 // index of next unused vertex
5786 unsigned int next_unused_vertex = 0;
5787
5788 // first the refinement of lines. children are stored
5789 // pairwise
5790 {
5791 // only active objects can be refined further
5794 endl = triangulation.end_line();
5796 next_unused_line = triangulation.begin_raw_line();
5797
5798 for (; line != endl; ++line)
5799 if (line->user_flag_set())
5800 {
5801 // this line needs to be refined
5802
5803 // find the next unused vertex and set it
5804 // appropriately
5805 while (triangulation.vertices_used[next_unused_vertex] == true)
5806 ++next_unused_vertex;
5807 Assert(
5808 next_unused_vertex < triangulation.vertices.size(),
5809 ExcMessage(
5810 "Internal error: During refinement, the triangulation wants to access an element of the 'vertices' array but it turns out that the array is not large enough."));
5811 triangulation.vertices_used[next_unused_vertex] = true;
5812
5813 triangulation.vertices[next_unused_vertex] = line->center(true);
5814
5815 // now that we created the right point, make up the
5816 // two child lines. To this end, find a pair of
5817 // unused lines
5818 bool pair_found = false;
5819 (void)pair_found;
5820 for (; next_unused_line != endl; ++next_unused_line)
5821 if (!next_unused_line->used() &&
5822 !(++next_unused_line)->used())
5823 {
5824 // go back to the first of the two unused
5825 // lines
5826 --next_unused_line;
5827 pair_found = true;
5828 break;
5829 }
5830 Assert(pair_found, ExcInternalError());
5831
5832 // there are now two consecutive unused lines, such
5833 // that the children of a line will be consecutive.
5834 // then set the child pointer of the present line
5835 line->set_children(0, next_unused_line->index());
5836
5837 // set the two new lines
5839 children[2] = {next_unused_line, ++next_unused_line};
5840 // some tests; if any of the iterators should be
5841 // invalid, then already dereferencing will fail
5842 AssertIsNotUsed(children[0]);
5843 AssertIsNotUsed(children[1]);
5844
5845 children[0]->set_bounding_object_indices(
5846 {line->vertex_index(0), next_unused_vertex});
5847 children[1]->set_bounding_object_indices(
5848 {next_unused_vertex, line->vertex_index(1)});
5849
5850 children[0]->set_used_flag();
5851 children[1]->set_used_flag();
5852 children[0]->clear_children();
5853 children[1]->clear_children();
5854 children[0]->clear_user_data();
5855 children[1]->clear_user_data();
5856 children[0]->clear_user_flag();
5857 children[1]->clear_user_flag();
5858
5859
5860 children[0]->set_boundary_id_internal(line->boundary_id());
5861 children[1]->set_boundary_id_internal(line->boundary_id());
5862
5863 children[0]->set_manifold_id(line->manifold_id());
5864 children[1]->set_manifold_id(line->manifold_id());
5865
5866 // finally clear flag indicating the need for
5867 // refinement
5868 line->clear_user_flag();
5869 }
5870 }
5871
5872
5873 // Now set up the new cells
5874
5875 // reserve space for inner lines (can be stored as single
5876 // lines)
5877 reserve_space(triangulation.faces->lines, 0, n_single_lines);
5878
5880 cells_with_distorted_children;
5881
5882 // reset next_unused_line, as now also single empty places in
5883 // the vector can be used
5885 next_unused_line = triangulation.begin_raw_line();
5886
5887 for (int level = 0;
5888 level < static_cast<int>(triangulation.levels.size()) - 1;
5889 ++level)
5890 {
5892 next_unused_cell = triangulation.begin_raw(level + 1);
5893
5894 for (const auto &cell :
5895 triangulation.active_cell_iterators_on_level(level))
5896 if (cell->refine_flag_set())
5897 {
5898 // actually set up the children and update neighbor
5899 // information
5900 create_children(triangulation,
5901 next_unused_vertex,
5902 next_unused_line,
5903 next_unused_cell,
5904 cell);
5905
5906 if (check_for_distorted_cells &&
5907 has_distorted_children<dim, spacedim>(cell))
5908 cells_with_distorted_children.distorted_cells.push_back(
5909 cell);
5910 // inform all listeners that cell refinement is done
5911 triangulation.signals.post_refinement_on_cell(cell);
5912 }
5913 }
5914
5915 return cells_with_distorted_children;
5916 }
5917
5918
5919 template <int spacedim>
5922 const bool check_for_distorted_cells)
5923 {
5924 static const int dim = 3;
5925 static const unsigned int X = numbers::invalid_unsigned_int;
5926 using raw_line_iterator =
5928 using raw_quad_iterator =
5930
5931 Assert(spacedim == 3, ExcNotImplemented());
5932
5933 Assert(triangulation.vertices.size() ==
5934 triangulation.vertices_used.size(),
5936
5937 // Check whether a new level is needed. We have to check for
5938 // this on the highest level only
5939 for (const auto &cell : triangulation.active_cell_iterators_on_level(
5940 triangulation.levels.size() - 1))
5941 if (cell->refine_flag_set())
5942 {
5943 triangulation.levels.push_back(
5944 std::make_unique<
5946 break;
5947 }
5948
5949 // first clear user flags for quads and lines; we're going to
5950 // use them to flag which lines and quads need refinement
5951 triangulation.faces->quads.clear_user_data();
5952 triangulation.faces->lines.clear_user_flags();
5953 triangulation.faces->quads.clear_user_flags();
5954
5955 // check how much space is needed on every level. We need not
5956 // check the highest level since either
5957 // - on the highest level no cells are flagged for refinement
5958 // - there are, but prepare_refinement added another empty
5959 // level which then is the highest level
5960
5961 // variables to hold the number of newly to be created
5962 // vertices, lines and quads. as these are stored globally,
5963 // declare them outside the loop over al levels. we need lines
5964 // and quads in pairs for refinement of old ones and lines and
5965 // quads, that can be stored as single ones, as they are newly
5966 // created in the inside of an existing cell
5967 unsigned int needed_vertices = 0;
5968 unsigned int needed_lines_single = 0;
5969 unsigned int needed_quads_single = 0;
5970 unsigned int needed_lines_pair = 0;
5971 unsigned int needed_quads_pair = 0;
5972 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
5973 {
5974 unsigned int new_cells = 0;
5975
5976 for (const auto &cell :
5977 triangulation.active_cell_iterators_on_level(level))
5978 if (cell->refine_flag_set())
5979 {
5980 // Only support isotropic refinement
5981 Assert(cell->refine_flag_set() ==
5984
5985 // Now count up how many new cells, faces, edges, and vertices
5986 // we will need to allocate to do this refinement.
5987 new_cells += cell->reference_cell().n_isotropic_children();
5988
5989 if (cell->reference_cell() == ReferenceCells::Hexahedron)
5990 {
5991 ++needed_vertices;
5992 needed_lines_single += 6;
5993 needed_quads_single += 12;
5994 }
5995 else if (cell->reference_cell() ==
5997 {
5998 needed_lines_single += 1;
5999 needed_quads_single += 8;
6000 }
6001 else
6002 {
6004 }
6005
6006 // Also check whether we have to refine any of the faces and
6007 // edges that bound this cell. They may of course already be
6008 // refined, so we only *mark* them for refinement by setting
6009 // the user flags
6010 for (const auto face : cell->face_indices())
6011 if (cell->face(face)->n_children() == 0)
6012 cell->face(face)->set_user_flag();
6013 else
6014 Assert(cell->face(face)->n_children() ==
6015 cell->reference_cell()
6016 .face_reference_cell(face)
6017 .n_isotropic_children(),
6019
6020 for (const auto line : cell->line_indices())
6021 if (cell->line(line)->has_children() == false)
6022 cell->line(line)->set_user_flag();
6023 else
6024 Assert(cell->line(line)->n_children() == 2,
6026 }
6027
6028 const unsigned int used_cells =
6029 std::count(triangulation.levels[level + 1]->cells.used.begin(),
6030 triangulation.levels[level + 1]->cells.used.end(),
6031 true);
6032
6033 if (triangulation.all_reference_cells_are_hyper_cube())
6034 reserve_space(*triangulation.levels[level + 1],
6035 used_cells + new_cells,
6036 3,
6037 spacedim,
6038 false);
6039 else
6040 reserve_space(*triangulation.levels[level + 1],
6041 used_cells + new_cells,
6042 3,
6043 spacedim,
6044 true);
6045
6046 reserve_space(triangulation.levels[level + 1]->cells, new_cells);
6047 }
6048
6049 // now count the quads and lines which were flagged for
6050 // refinement
6053 quad != triangulation.end_quad();
6054 ++quad)
6055 {
6056 if (quad->user_flag_set() == false)
6057 continue;
6058
6059 if (quad->reference_cell() == ReferenceCells::Quadrilateral)
6060 {
6061 needed_quads_pair += 4;
6062 needed_lines_pair += 4;
6063 needed_vertices += 1;
6064 }
6065 else if (quad->reference_cell() == ReferenceCells::Triangle)
6066 {
6067 needed_quads_pair += 4;
6068 needed_lines_single += 3;
6069 }
6070 else
6071 {
6073 }
6074 }
6075
6078 line != triangulation.end_line();
6079 ++line)
6080 {
6081 if (line->user_flag_set() == false)
6082 continue;
6083
6084 needed_lines_pair += 2;
6085 needed_vertices += 1;
6086 }
6087
6088 reserve_space(triangulation.faces->lines,
6089 needed_lines_pair,
6090 needed_lines_single);
6092 needed_quads_pair,
6093 needed_quads_single);
6094 reserve_space(triangulation.faces->quads,
6095 needed_quads_pair,
6096 needed_quads_single);
6097
6098
6099 // add to needed vertices how many vertices are already in use
6100 needed_vertices += std::count(triangulation.vertices_used.begin(),
6101 triangulation.vertices_used.end(),
6102 true);
6103
6104 if (needed_vertices > triangulation.vertices.size())
6105 {
6106 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
6107 triangulation.vertices_used.resize(needed_vertices, false);
6108 }
6109
6110 //-----------------------------------------
6111 // Before we start with the actual refinement, we do some
6112 // sanity checks if in debug mode. especially, we try to catch
6113 // the notorious problem with lines being twice refined,
6114 // i.e. there are cells adjacent at one line ("around the
6115 // edge", but not at a face), with two cells differing by more
6116 // than one refinement level
6117 //
6118 // this check is very simple to implement here, since we have
6119 // all lines flagged if they shall be refined
6120#ifdef DEBUG
6121 for (const auto &cell : triangulation.active_cell_iterators())
6122 if (!cell->refine_flag_set())
6123 for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
6124 if (cell->line(line_n)->has_children())
6125 for (unsigned int c = 0; c < 2; ++c)
6126 Assert(cell->line(line_n)->child(c)->user_flag_set() == false,
6128#endif
6129
6130 unsigned int current_vertex = 0;
6131
6132 // helper function - find the next available vertex number and mark it
6133 // as used.
6134 auto get_next_unused_vertex = [](const unsigned int current_vertex,
6135 std::vector<bool> &vertices_used) {
6136 unsigned int next_vertex = current_vertex;
6137 while (next_vertex < vertices_used.size() &&
6138 vertices_used[next_vertex] == true)
6139 ++next_vertex;
6140 Assert(next_vertex < vertices_used.size(), ExcInternalError());
6141 vertices_used[next_vertex] = true;
6142
6143 return next_vertex;
6144 };
6145
6146 // LINES
6147 {
6150 endl = triangulation.end_line();
6151 raw_line_iterator next_unused_line = triangulation.begin_raw_line();
6152
6153 for (; line != endl; ++line)
6154 {
6155 if (line->user_flag_set() == false)
6156 continue;
6157
6158 next_unused_line =
6159 triangulation.faces->lines.template next_free_pair_object<1>(
6161 Assert(next_unused_line.state() == IteratorState::valid,
6163
6164 // now we found two consecutive unused lines, such
6165 // that the children of a line will be consecutive.
6166 // then set the child pointer of the present line
6167 line->set_children(0, next_unused_line->index());
6168
6169 const std::array<raw_line_iterator, 2> children{
6170 {next_unused_line, ++next_unused_line}};
6171
6172 AssertIsNotUsed(children[0]);
6173 AssertIsNotUsed(children[1]);
6174
6175 current_vertex =
6176 get_next_unused_vertex(current_vertex,
6177 triangulation.vertices_used);
6178 triangulation.vertices[current_vertex] = line->center(true);
6179
6180 children[0]->set_bounding_object_indices(
6181 {line->vertex_index(0), current_vertex});
6182 children[1]->set_bounding_object_indices(
6183 {current_vertex, line->vertex_index(1)});
6184
6185 const auto manifold_id = line->manifold_id();
6186 const auto boundary_id = line->boundary_id();
6187 for (const auto &child : children)
6188 {
6189 child->set_used_flag();
6190 child->clear_children();
6191 child->clear_user_data();
6192 child->clear_user_flag();
6193 child->set_boundary_id_internal(boundary_id);
6194 child->set_manifold_id(manifold_id);
6195 }
6196
6197 line->clear_user_flag();
6198 }
6199 }
6200
6201 // QUADS
6202 {
6204 quad = triangulation.begin_quad(),
6205 endq = triangulation.end_quad();
6206
6207 for (; quad != endq; ++quad)
6208 {
6209 if (quad->user_flag_set() == false)
6210 continue;
6211
6212 const auto reference_face_type = quad->reference_cell();
6213
6214 // 1) create new lines (property is set later)
6215 // maximum of 4 new lines (4 quadrilateral, 3 triangle)
6216 std::array<raw_line_iterator, 4> new_lines;
6217 if (reference_face_type == ReferenceCells::Quadrilateral)
6218 {
6219 for (unsigned int l = 0; l < 2; ++l)
6220 {