From 9aa53960a5dbdfb9ea2418916871d8c403b6feb2 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Fri, 20 Sep 2024 13:58:50 +0200 Subject: [PATCH] Add support for simplex meshes (#132) --- examples/benchmarks_3D.cc | 34 +++-- examples/repairing.cc | 232 ++++++++++++++++++++++++++++++++ include/agglomeration_handler.h | 20 ++- include/poly_utils.h | 44 ++++-- source/agglomeration_handler.cc | 3 +- 5 files changed, 303 insertions(+), 30 deletions(-) create mode 100644 examples/repairing.cc diff --git a/examples/benchmarks_3D.cc b/examples/benchmarks_3D.cc index 4ae08298..c9636797 100644 --- a/examples/benchmarks_3D.cc +++ b/examples/benchmarks_3D.cc @@ -1,3 +1,6 @@ +#include +#include + #include #include #include @@ -47,7 +50,7 @@ class AgglomerationBenchmark make_grid(); Triangulation tria; - MappingQ mapping; + MappingFE mapping; std::unique_ptr> ah; public: @@ -74,7 +77,7 @@ AgglomerationBenchmark::AgglomerationBenchmark( const PartitionerType &partitioner_type, const unsigned int extraction_level, const unsigned int n_subdomains) - : mapping(1) + : mapping(FE_SimplexDGP(1)) , grid_type(grid_type) , partitioner_type(partitioner_type) , extraction_level(extraction_level) @@ -105,9 +108,11 @@ AgglomerationBenchmark::make_grid() // } grid_in.attach_triangulation(tria); std::ifstream mesh_file( - "../../meshes/csf_brain_filled_centered_UCD.inp"); // piston mesh - grid_in.read_ucd(mesh_file); - tria.refine_global(1); + "../../meshes/gray_level_image1.vtk"); // liver mesh + grid_in.read_vtk(mesh_file); + // grid_in.read_ucd(mesh_file); + // "../../meshes/csf_brain_filled_centered_UCD.inp"); // brain mesh + // tria.refine_global(1); } } else @@ -158,8 +163,8 @@ AgglomerationBenchmark::make_grid() // const auto tree = pack_rtree>(boxes); - auto start = std::chrono::system_clock::now(); - const auto tree = pack_rtree>(boxes); + auto start = std::chrono::system_clock::now(); + auto tree = pack_rtree>(boxes); std::cout << "Total number of available levels: " << n_levels(tree) << std::endl; @@ -169,9 +174,9 @@ AgglomerationBenchmark::make_grid() ExcMessage("At least two levels are needed.")); #endif - const auto &csr_and_agglomerates = - PolyUtils::extract_children_of_level(tree, extraction_level); - const auto &agglomerates = csr_and_agglomerates.second; // vec> + CellsAgglomerator agglomerator{tree, + extraction_level}; + const auto agglomerates = agglomerator.extract_agglomerates(); std::size_t agglo_index = 0; for (std::size_t i = 0; i < agglomerates.size(); ++i) @@ -240,7 +245,7 @@ main() std::cout << "Benchmarking with Rtree:" << std::endl; // for (unsigned int i = 0; i < 10; ++i) // { - for (const unsigned int extraction_level : {2, 3, 4, 5}) + for (const unsigned int extraction_level : {2, 3, 4, 5, 6}) { AgglomerationBenchmark<3> poisson_problem{GridType::unstructured, PartitionerType::rtree, @@ -255,12 +260,15 @@ main() // // { // // for (const unsigned int target_partitions : {12, 90, 715, 5715}) // //piston - for (const unsigned int target_partitions : {47, 372, 2976, 23804}) // brain + // for (const unsigned int target_partitions : {47, 372, 2976, 23804}) // + // brain + for (const unsigned int target_partitions : + {9, 70, 556, 4441, 29000}) // liver { AgglomerationBenchmark<3> poisson_problem{GridType::unstructured, PartitionerType::metis, 0, - target_partitions}; // 16, 64 + target_partitions}; poisson_problem.run(); } } diff --git a/examples/repairing.cc b/examples/repairing.cc new file mode 100644 index 00000000..e8517a95 --- /dev/null +++ b/examples/repairing.cc @@ -0,0 +1,232 @@ +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include + + +// Check that the R-tree based agglomeration works also in 3D. + + +enum class GridType +{ + grid_generator, + unstructured +}; + +enum class PartitionerType +{ + metis, + rtree +}; + + +template +class Test +{ +private: + void + make_grid(); + void + setup_agglomeration(); + + Triangulation tria; + MappingFE mapping; + FE_DGQ dg_fe; + std::unique_ptr> ah; + // no hanging node in DG discretization, we define an AffineConstraints object + // so we can use the distribute_local_to_global() directly. + AffineConstraints constraints; + Vector system_rhs; + std::unique_ptr> cached_tria; + std::unique_ptr> rhs_function; + +public: + Test(const GridType &grid_type = GridType::grid_generator, + const PartitionerType &partitioner_type = PartitionerType::rtree, + const unsigned int = 0, + const unsigned int = 0, + const unsigned int fe_degree = 1); + void + run(); + + double penalty_constant = 10.; + GridType grid_type; + PartitionerType partitioner_type; + unsigned int extraction_level; + unsigned int n_subdomains; +}; + + + +template +Test::Test(const GridType &grid_type, + const PartitionerType &partitioner_type, + const unsigned int extraction_level, + const unsigned int n_subdomains, + const unsigned int fe_degree) + : mapping(FE_SimplexDGP(fe_degree)) + , dg_fe(fe_degree) + , grid_type(grid_type) + , partitioner_type(partitioner_type) + , extraction_level(extraction_level) + , n_subdomains(n_subdomains) +{} + +template +void +Test::make_grid() +{ + GridIn grid_in; + if (grid_type == GridType::unstructured) + { + grid_in.attach_triangulation(tria); + std::cout << "########### Reading mesh file... ###########" << std::endl; + std::ifstream filename("../../meshes/ernie.msh"); // liver or brain domain + grid_in.read_msh(filename); + std::cout << "########### Done ###########" << std::endl; + } + else + { + GridGenerator::hyper_cube(tria, 0, 1); + tria.refine_global(2); + GridTools::distort_random(0.25, tria); + } + std::cout << "Size of tria: " << tria.n_active_cells() << std::endl; + cached_tria = std::make_unique>(tria, mapping); +} + +template +void +Test::setup_agglomeration() + +{ + ah = std::make_unique>(*cached_tria); + + std::vector idxs_to_be_agglomerated = {0, 1, 2, 7}; + + std::vector::active_cell_iterator> + cells_to_be_agglomerated; + PolyUtils::collect_cells_for_agglomeration(tria, + idxs_to_be_agglomerated, + cells_to_be_agglomerated); + + std::vector idxs_to_be_agglomerated2 = {4, 5, 9}; + + std::vector::active_cell_iterator> + cells_to_be_agglomerated2; + PolyUtils::collect_cells_for_agglomeration(tria, + idxs_to_be_agglomerated2, + cells_to_be_agglomerated2); + + std::vector idxs_to_be_agglomerated3 = {3, 6}; + + std::vector::active_cell_iterator> + cells_to_be_agglomerated3; + PolyUtils::collect_cells_for_agglomeration(tria, + idxs_to_be_agglomerated3, + cells_to_be_agglomerated3); + + std::vector idxs_to_be_agglomerated4 = {8, + 9, + 10, + 15}; + + std::vector::active_cell_iterator> + cells_to_be_agglomerated4; + PolyUtils::collect_cells_for_agglomeration(tria, + idxs_to_be_agglomerated4, + cells_to_be_agglomerated4); + + std::vector idxs_to_be_agglomerated5 = {11, + 12, + 13, + 14}; + + std::vector::active_cell_iterator> + cells_to_be_agglomerated5; + PolyUtils::collect_cells_for_agglomeration(tria, + idxs_to_be_agglomerated5, + cells_to_be_agglomerated5); + + // Agglomerate the cells just stored + ah->define_agglomerate_with_check(cells_to_be_agglomerated); + ah->define_agglomerate_with_check(cells_to_be_agglomerated2); + ah->define_agglomerate_with_check(cells_to_be_agglomerated3); + ah->define_agglomerate_with_check(cells_to_be_agglomerated4); + ah->define_agglomerate_with_check(cells_to_be_agglomerated5); + + std::cout << "Number of generated agglomerates: " << ah->n_agglomerates() + << std::endl; + + // std::cout << "########### Repairing the grid... ###########" << + // std::endl; ah->repair_grid(); + std::cout << "Defined all agglomerates" << std::endl; + + // Check local bboxes + for (const auto &box : ah->get_local_bboxes()) + { + std::cout << "p0: " << box.get_boundary_points().first << std::endl; + std::cout << "p1: " << box.get_boundary_points().second << std::endl; + std::cout << std::endl; + } + + std::cout << "########### Create output for visualization... ###########" + << std::endl; + { + const std::string &partitioner = + (partitioner_type == PartitionerType::metis) ? "metis" : "rtree"; + + const std::string filename = "square_repaired" + partitioner + ".vtu"; + std::ofstream output(filename); + + DataOut data_out; + data_out.attach_dof_handler(ah->agglo_dh); + // data_out.attach_dof_handler(ah->output_dh); + + Vector agglo_idx(tria.n_active_cells()); + for (const auto &polytope : ah->polytope_iterators()) + { + const types::global_cell_index polytopal_idx = polytope->index(); + const auto &deal_cells = polytope->get_agglomerate(); + for (const auto &cell : deal_cells) + agglo_idx[cell->active_cell_index()] = polytopal_idx; + } + data_out.add_data_vector(agglo_idx, + "agglomerated_idx", + DataOut::type_cell_data); + data_out.build_patches(mapping); + data_out.write_vtu(output); + } + std::cout << "########### Done ###########" << std::endl; +} + +template +void +Test::run() +{ + make_grid(); + setup_agglomeration(); +} + +int +main() +{ + { + Test<2> test{GridType::grid_generator, + PartitionerType::rtree, + 4 /*extraction_level*/}; + test.run(); + } + + + + return 0; +} diff --git a/include/agglomeration_handler.h b/include/agglomeration_handler.h index b273e6d3..c3c0748a 100644 --- a/include/agglomeration_handler.h +++ b/include/agglomeration_handler.h @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -56,9 +57,6 @@ using namespace dealii; // Forward declarations -template -class FinitElement; - template class AgglomerationHandler; @@ -298,6 +296,9 @@ class AgglomerationHandler : public Subscriptor inline const FiniteElement & get_fe() const; + inline const Mapping & + get_mapping() const; + inline const std::vector> & get_local_bboxes() const; @@ -831,6 +832,10 @@ class AgglomerationHandler : public Subscriptor // Map the master cell index with the polytope index std::map master2polygon; + + std::vector::active_cell_iterator> + master_disconnected; + // Dummy FiniteElement objects needed only to generate quadratures /** @@ -883,6 +888,15 @@ AgglomerationHandler::get_fe() const +template +inline const Mapping & +AgglomerationHandler::get_mapping() const +{ + return *mapping; +} + + + template inline const Triangulation & AgglomerationHandler::get_triangulation() const diff --git a/include/poly_utils.h b/include/poly_utils.h index ae9a7f7e..f14c5869 100644 --- a/include/poly_utils.h +++ b/include/poly_utils.h @@ -878,13 +878,22 @@ namespace dealii::PolyUtils const_cast *>( &agglomeration_handler.output_dh); const FiniteElement &fe = agglomeration_handler.get_fe(); + const Mapping &mapping = agglomeration_handler.get_mapping(); const Triangulation &tria = agglomeration_handler.get_triangulation(); const auto &bboxes = agglomeration_handler.get_local_bboxes(); + std::unique_ptr> output_fe; + if (tria.all_reference_cells_are_hyper_cube()) + output_fe = std::make_unique>(fe.degree); + else if (tria.all_reference_cells_are_simplex()) + output_fe = std::make_unique>(fe.degree); + else + AssertThrow(false, ExcNotImplemented()); + // Setup an auxiliary DoFHandler for output purposes output_dh->reinit(tria); - output_dh->distribute_dofs(fe); + output_dh->distribute_dofs(*output_fe); const IndexSet &locally_owned_dofs = output_dh->locally_owned_dofs(); const IndexSet locally_relevant_dofs = @@ -900,10 +909,12 @@ namespace dealii::PolyUtils std::vector agglo_dof_indices(fe.dofs_per_cell); std::vector standard_dof_indices( fe.dofs_per_cell); - std::vector output_dof_indices(fe.dofs_per_cell); + std::vector output_dof_indices( + output_fe->dofs_per_cell); - Quadrature quad(fe.get_unit_support_points()); - FEValues output_fe_values(fe, + Quadrature quad(output_fe->get_unit_support_points()); + FEValues output_fe_values(mapping, + *output_fe, quad, update_quadrature_points); @@ -1076,19 +1087,25 @@ namespace dealii::PolyUtils // otherwise, do not create any matrix const Triangulation &tria = agglomeration_handler.get_triangulation(); + const Mapping &mapping = agglomeration_handler.get_mapping(); const FiniteElement &original_fe = agglomeration_handler.get_fe(); - // We use DGQ nodal elements of the same degree as the ones in the - // agglomeration handler to generate the output also in the case in - // which different elements are used, such as DGP. + // We use DGQ (on tensor-product meshes) or DGP (on simplex meshes) + // nodal elements of the same degree as the ones in the agglomeration + // handler to interpolate the solution onto the finer grid. + std::unique_ptr> output_fe; + if (tria.all_reference_cells_are_hyper_cube()) + output_fe = std::make_unique>(original_fe.degree); + else if (tria.all_reference_cells_are_simplex()) + output_fe = std::make_unique>(original_fe.degree); + else + AssertThrow(false, ExcNotImplemented()); - FE_DGQ output_fe(original_fe.degree); DoFHandler &output_dh = const_cast &>(agglomeration_handler.output_dh); - output_dh.reinit(tria); - output_dh.distribute_dofs(output_fe); + output_dh.distribute_dofs(*output_fe); if constexpr (std::is_same_v) { @@ -1112,9 +1129,10 @@ namespace dealii::PolyUtils const unsigned int dofs_per_cell = agglomeration_handler.n_dofs_per_cell(); - const unsigned int output_dofs_per_cell = output_fe.n_dofs_per_cell(); - Quadrature quad(output_fe.get_unit_support_points()); - FEValues output_fe_values(output_fe, + const unsigned int output_dofs_per_cell = output_fe->n_dofs_per_cell(); + Quadrature quad(output_fe->get_unit_support_points()); + FEValues output_fe_values(mapping, + *output_fe, quad, update_quadrature_points); diff --git a/source/agglomeration_handler.cc b/source/agglomeration_handler.cc index 5eda642c..a412d103 100644 --- a/source/agglomeration_handler.cc +++ b/source/agglomeration_handler.cc @@ -142,7 +142,6 @@ AgglomerationHandler::define_agglomerate_with_check( } - template void AgglomerationHandler::initialize_fe_values( @@ -247,6 +246,8 @@ AgglomerationHandler::distribute_agglomerated_dofs( fe = std::make_unique>(fe_space.degree); else if (dynamic_cast *>(&fe_space)) fe = std::make_unique>(fe_space.degree); + else if (dynamic_cast *>(&fe_space)) + fe = std::make_unique>(fe_space.degree); else AssertThrow( false,