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