22#ifndef INCLUDE_EXADG_GRID_GRID_UTILITIES_H_
23#define INCLUDE_EXADG_GRID_GRID_UTILITIES_H_
26#include <deal.II/fe/fe_simplex_p.h>
27#include <deal.II/fe/mapping_fe.h>
28#include <deal.II/fe/mapping_q.h>
29#include <deal.II/grid/grid_in.h>
30#include <deal.II/grid/grid_tools.h>
31#include <deal.II/multigrid/mg_transfer_global_coarsening.h>
34#include <exadg/grid/balanced_granularity_partition_policy.h>
35#include <exadg/grid/grid.h>
36#include <exadg/grid/grid_data.h>
37#include <exadg/grid/perform_local_refinements.h>
41namespace GridUtilities
44using PeriodicFacePairs = std::vector<
45 dealii::GridTools::PeriodicFacePair<typename dealii::Triangulation<dim>::cell_iterator>>;
52create_mapping(std::shared_ptr<dealii::Mapping<dim>> & mapping,
53 ElementType
const & element_type,
54 unsigned int const & mapping_degree)
56 if(element_type == ElementType::Hypercube)
58 mapping = std::make_shared<dealii::MappingQ<dim>>(mapping_degree);
60 else if(element_type == ElementType::Simplex)
62 mapping = std::make_shared<dealii::MappingFE<dim>>(dealii::FE_SimplexP<dim>(mapping_degree));
66 AssertThrow(
false, dealii::ExcMessage(
"Invalid parameter element_type."));
76template<
int dim,
typename Number>
78create_mapping_with_multigrid(std::shared_ptr<dealii::Mapping<dim>> & mapping,
79 std::shared_ptr<MultigridMappings<dim, Number>> & multigrid_mappings,
80 ElementType
const & element_type,
81 unsigned int const & mapping_degree_fine,
82 unsigned int const & mapping_degree_coarse,
83 bool const involves_h_multigrid)
86 create_mapping(mapping, element_type, mapping_degree_fine);
89 std::shared_ptr<dealii::Mapping<dim>> coarse_mapping;
90 if(involves_h_multigrid and (mapping_degree_coarse != mapping_degree_fine))
92 create_mapping(coarse_mapping, element_type, mapping_degree_coarse);
94 multigrid_mappings = std::make_shared<MultigridMappings<dim, Number>>(mapping, coarse_mapping);
103std::vector<dealii::GridTools::PeriodicFacePair<typename dealii::DoFHandler<dim>::cell_iterator>>
104transform_periodic_face_pairs_to_dof_cell_iterator(
105 std::vector<dealii::GridTools::PeriodicFacePair<
106 typename dealii::Triangulation<dim>::cell_iterator>>
const & periodic_faces,
107 dealii::DoFHandler<dim>
const & dof_handler)
109 typedef typename std::vector<
110 dealii::GridTools::PeriodicFacePair<typename dealii::DoFHandler<dim>::cell_iterator>>
113 PeriodicFacesDoF periodic_faces_dof;
115 for(
auto it : periodic_faces)
117 dealii::GridTools::PeriodicFacePair<typename dealii::DoFHandler<dim>::cell_iterator>
118 face_pair_dof_hander;
120 face_pair_dof_hander.cell[0] = it.cell[0]->as_dof_handler_iterator(dof_handler);
121 face_pair_dof_hander.cell[1] = it.cell[1]->as_dof_handler_iterator(dof_handler);
123 face_pair_dof_hander.face_idx[0] = it.face_idx[0];
124 face_pair_dof_hander.face_idx[1] = it.face_idx[1];
126 face_pair_dof_hander.orientation = it.orientation;
127 face_pair_dof_hander.matrix = it.matrix;
129 periodic_faces_dof.push_back(face_pair_dof_hander);
132 return periodic_faces_dof;
146 std::shared_ptr<dealii::Triangulation<dim>> & triangulation,
147 PeriodicFacePairs<dim> & periodic_face_pairs,
148 MPI_Comm
const & mpi_comm,
149 GridData
const & data,
150 bool const construct_multigrid_hierarchy,
151 std::function<
void(dealii::Triangulation<dim> &,
152 PeriodicFacePairs<dim> &,
154 std::vector<unsigned int>
const &)>
const & lambda_create_triangulation,
155 unsigned int const global_refinements,
156 std::vector<unsigned int>
const & vector_local_refinements)
158 if(vector_local_refinements.size() != 0)
161 data.element_type == ElementType::Hypercube,
163 "Local refinements are currently only supported for meshes composed of hypercube elements."));
167 auto mesh_smoothing = dealii::Triangulation<dim>::none;
173 if(data.element_type == ElementType::Hypercube)
175 mesh_smoothing = dealii::Triangulation<dim>::limit_level_difference_at_vertices;
178 if(data.triangulation_type == TriangulationType::Serial)
180 AssertDimension(dealii::Utilities::MPI::n_mpi_processes(mpi_comm), 1);
181 triangulation = std::make_shared<dealii::Triangulation<dim>>(mesh_smoothing);
183 lambda_create_triangulation(*triangulation,
186 vector_local_refinements);
188 else if(data.triangulation_type == TriangulationType::Distributed)
190 typename dealii::parallel::distributed::Triangulation<dim>::Settings distributed_settings;
192 if(construct_multigrid_hierarchy)
194 distributed_settings =
195 dealii::parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy;
199 std::make_shared<dealii::parallel::distributed::Triangulation<dim>>(mpi_comm,
201 distributed_settings);
203 lambda_create_triangulation(*triangulation,
206 vector_local_refinements);
208 else if(data.triangulation_type == TriangulationType::FullyDistributed)
210 auto const serial_grid_generator = [&](dealii::Triangulation<dim, dim> & tria_serial) {
211 lambda_create_triangulation(tria_serial,
214 vector_local_refinements);
217 auto const serial_grid_partitioner = [&](dealii::Triangulation<dim, dim> & tria_serial,
219 unsigned int const group_size) {
221 if(data.partitioning_type == PartitioningType::Metis)
223 dealii::GridTools::partition_triangulation(dealii::Utilities::MPI::n_mpi_processes(comm),
226 else if(data.partitioning_type == PartitioningType::z_order)
228 dealii::GridTools::partition_triangulation_zorder(
229 dealii::Utilities::MPI::n_mpi_processes(comm), tria_serial);
233 AssertThrow(
false, dealii::ExcNotImplemented());
237 unsigned int const group_size = 1;
239 typename dealii::TriangulationDescription::Settings triangulation_description_setting =
240 dealii::TriangulationDescription::default_setting;
242 if(construct_multigrid_hierarchy)
244 triangulation_description_setting =
245 dealii::TriangulationDescription::construct_multigrid_hierarchy;
249 std::make_shared<dealii::parallel::fullydistributed::Triangulation<dim>>(mpi_comm);
251 auto const description = dealii::TriangulationDescription::Utilities::
252 create_description_from_triangulation_in_groups<dim, dim>(serial_grid_generator,
253 serial_grid_partitioner,
254 triangulation->get_communicator(),
257 triangulation_description_setting);
259 triangulation->create_triangulation(description);
263 AssertThrow(
false, dealii::ExcMessage(
"Invalid parameter triangulation_type."));
280create_coarse_triangulations_automatically_from_fine_triangulation(
281 dealii::Triangulation<dim>
const & fine_triangulation,
282 PeriodicFacePairs<dim>
const & fine_periodic_face_pairs,
283 std::vector<std::shared_ptr<dealii::Triangulation<dim>
const>> & coarse_triangulations_const,
284 std::vector<PeriodicFacePairs<dim>> & coarse_periodic_face_pairs,
285 GridData
const & data)
289 AssertThrow(data.triangulation_type == TriangulationType::Serial or
290 data.triangulation_type == TriangulationType::Distributed,
291 dealii::ExcMessage(
"Invalid parameter triangulation_type."));
293 if(data.triangulation_type == TriangulationType::Serial)
296 fine_triangulation.all_reference_cells_are_hyper_cube(),
298 "The create_geometric_coarsening_sequence function of dealii does currently not support "
299 "simplicial elements."));
301 coarse_triangulations_const =
302 dealii::MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence(
305 else if(data.triangulation_type == TriangulationType::Distributed)
308 fine_triangulation.all_reference_cells_are_hyper_cube(),
310 "dealii::parallel::distributed::Triangulation does not support simplicial elements."));
312 coarse_triangulations_const =
313 dealii::MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence(
315 BalancedGranularityPartitionPolicy<dim>(
316 dealii::Utilities::MPI::n_mpi_processes(fine_triangulation.get_communicator())));
320 AssertThrow(
false, dealii::ExcMessage(
"Invalid parameter triangulation_type."));
326 coarse_triangulations_const.pop_back();
328 coarse_periodic_face_pairs.resize(coarse_triangulations_const.size());
329 for(
unsigned int level = 0; level < coarse_periodic_face_pairs.size(); ++level)
331 coarse_periodic_face_pairs[level] = fine_periodic_face_pairs;
337create_coarse_triangulations_for_fully_distributed_triangulation(
338 dealii::Triangulation<dim>
const & fine_triangulation,
339 std::vector<std::shared_ptr<dealii::Triangulation<dim>
const>> & coarse_triangulations_const,
340 std::vector<PeriodicFacePairs<dim>> & coarse_periodic_face_pairs,
341 GridData
const & data,
342 std::function<
void(dealii::Triangulation<dim> &,
343 PeriodicFacePairs<dim> &,
345 std::vector<unsigned int>
const &)>
const & lambda_create_triangulation,
346 std::vector<unsigned int>
const vector_local_refinements)
351 AssertThrow(data.triangulation_type == TriangulationType::FullyDistributed,
352 dealii::ExcMessage(
"Invalid parameter triangulation_type."));
354 if(fine_triangulation.n_global_levels() >= 2)
357 std::vector<std::shared_ptr<dealii::Triangulation<dim>>> coarse_triangulations =
358 std::vector<std::shared_ptr<dealii::Triangulation<dim>>>(
359 fine_triangulation.n_global_levels() - 1);
361 coarse_periodic_face_pairs = std::vector<PeriodicFacePairs<dim>>(coarse_triangulations.size());
364 unsigned int level = fine_triangulation.n_global_levels() - 2;
365 std::vector<unsigned int>
refine_local = vector_local_refinements;
368 if(data.n_refine_global >= 1)
370 unsigned int const n_refine_global_start = (
unsigned int)(data.n_refine_global - 1);
371 for(
int refine_global = n_refine_global_start; refine_global >= 0; --refine_global)
373 GridUtilities::create_triangulation<dim>(coarse_triangulations[level],
374 coarse_periodic_face_pairs[level],
375 fine_triangulation.get_communicator(),
378 lambda_create_triangulation,
394 for(
size_t material_id = 0; material_id <
refine_local.size(); material_id++)
402 GridUtilities::create_triangulation<dim>(coarse_triangulations[level],
403 coarse_periodic_face_pairs[level],
404 fine_triangulation.get_communicator(),
407 lambda_create_triangulation,
421 "There occurred a logical error when creating the geometric coarsening sequence."));
424 for(
auto const & it : coarse_triangulations)
426 coarse_triangulations_const.push_back(it);
433create_coarse_triangulations(
434 dealii::Triangulation<dim>
const & fine_triangulation,
435 PeriodicFacePairs<dim>
const & fine_periodic_face_pairs,
436 std::vector<std::shared_ptr<dealii::Triangulation<dim>
const>> & coarse_triangulations_const,
437 std::vector<PeriodicFacePairs<dim>> & coarse_periodic_face_pairs,
438 GridData
const & data,
439 std::function<
void(dealii::Triangulation<dim> &,
440 PeriodicFacePairs<dim> &,
442 std::vector<unsigned int>
const &)>
const & lambda_create_triangulation,
443 std::vector<unsigned int>
const vector_local_refinements)
448 if(data.triangulation_type == TriangulationType::Serial or
449 data.triangulation_type == TriangulationType::Distributed)
451 create_coarse_triangulations_automatically_from_fine_triangulation(fine_triangulation,
452 fine_periodic_face_pairs,
453 coarse_triangulations_const,
454 coarse_periodic_face_pairs,
457 else if(data.triangulation_type == TriangulationType::FullyDistributed)
459 create_coarse_triangulations_for_fully_distributed_triangulation(fine_triangulation,
460 coarse_triangulations_const,
461 coarse_periodic_face_pairs,
463 lambda_create_triangulation,
464 vector_local_refinements);
468 AssertThrow(
false, dealii::ExcMessage(
"Invalid parameter triangulation_type."));
482 MPI_Comm
const & mpi_comm,
483 GridData
const & data,
484 std::function<
void(dealii::Triangulation<dim> &,
485 PeriodicFacePairs<dim> &,
487 std::vector<unsigned int>
const &)>
const & lambda_create_triangulation,
488 std::vector<unsigned int>
const vector_local_refinements)
490 GridUtilities::create_triangulation(grid.triangulation,
491 grid.periodic_face_pairs,
495 lambda_create_triangulation,
496 data.n_refine_global,
497 vector_local_refinements);
507create_triangulation_with_multigrid(
509 MPI_Comm
const & mpi_comm,
510 GridData
const & data,
511 bool const involves_h_multigrid,
512 std::function<
void(dealii::Triangulation<dim> &,
513 PeriodicFacePairs<dim> &,
515 std::vector<unsigned int>
const &)>
const & lambda_create_triangulation,
516 std::vector<unsigned int>
const vector_local_refinements)
518 if(involves_h_multigrid)
521 if(data.element_type == ElementType::Simplex)
523 AssertThrow(data.create_coarse_triangulations ==
true,
525 "You need to set GridData::create_coarse_triangulations = true "
526 "in order to use h-multigrid for simplex meshes."));
530 GridUtilities::create_triangulation(grid.triangulation,
531 grid.periodic_face_pairs,
534 not data.create_coarse_triangulations,
535 lambda_create_triangulation,
536 data.n_refine_global,
537 vector_local_refinements);
540 if(grid.triangulation->has_hanging_nodes())
542 AssertThrow(data.create_coarse_triangulations ==
true,
544 "You need to set GridData::create_coarse_triangulations = true "
545 "in order to use h-multigrid for meshes with hanging nodes."));
549 if(data.create_coarse_triangulations)
551 GridUtilities::create_coarse_triangulations(*grid.triangulation,
552 grid.periodic_face_pairs,
553 grid.coarse_triangulations,
554 grid.coarse_periodic_face_pairs,
556 lambda_create_triangulation,
557 vector_local_refinements);
564 GridUtilities::create_triangulation<dim>(
565 grid, mpi_comm, data, lambda_create_triangulation, vector_local_refinements);
575create_coarse_triangulations_after_coarsening_and_refinement(
576 dealii::Triangulation<dim>
const & fine_triangulation,
577 PeriodicFacePairs<dim>
const & fine_periodic_face_pairs,
578 std::vector<std::shared_ptr<dealii::Triangulation<dim>
const>> & coarse_triangulations_const,
579 std::vector<PeriodicFacePairs<dim>> & coarse_periodic_face_pairs,
580 GridData
const & data,
581 bool const amr_preserves_boundary_cells)
583 if(data.triangulation_type == TriangulationType::Serial or
584 data.triangulation_type == TriangulationType::Distributed)
587 if(fine_periodic_face_pairs.size() > 0)
589 AssertThrow(amr_preserves_boundary_cells,
591 "Combination of adaptive mesh refinement and periodic face pairs"
592 " requires boundary cells to be preserved."));
595 create_coarse_triangulations_automatically_from_fine_triangulation(fine_triangulation,
596 fine_periodic_face_pairs,
597 coarse_triangulations_const,
598 coarse_periodic_face_pairs,
601 else if(data.triangulation_type == TriangulationType::FullyDistributed)
604 dealii::ExcMessage(
"Combination of adaptive mesh refinement and "
605 "TriangulationType::FullyDistributed not implemented."));
609 AssertThrow(
false, dealii::ExcMessage(
"Invalid parameter triangulation_type."));
619read_external_triangulation(dealii::Triangulation<dim, dim> & tria, GridData
const & data)
621 AssertThrow(not data.file_name.empty(),
623 "You are trying to read a grid file, but the string, which is supposed to contain"
624 " the file name, is empty."));
626 dealii::GridIn<dim> grid_in;
628 grid_in.attach_triangulation(tria);
631 std::string extension = data.file_name.substr(data.file_name.find_last_of(
'.') + 1);
633 AssertThrow(not extension.empty(),
634 dealii::ExcMessage(
"You are trying to read a grid file, but the file extension is"
638 typename dealii::GridIn<dim>::Format format;
639 if(extension ==
"e" || extension ==
"exo")
640 format = dealii::GridIn<dim>::Format::exodusii;
642 format = grid_in.parse_format(extension);
647 grid_in.read(data.file_name, format);
650 dealii::ExcMessage(
"You are trying to read a grid file, but the element type of the"
651 " external grid file and the element type specified in GridData"
652 " don't match. Most likely, you forgot to change the element_type"
653 " parameter of GridData to the desired element type."));
ElementType get_element_type(dealii::Triangulation< dim > const &tria)
Definition grid_data.h:70
void refine_local(dealii::Triangulation< dim > &tria, std::vector< unsigned int > const &vector_local_refinements)
Definition perform_local_refinements.h:43