Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Enable simplex meshes #132

Merged
merged 1 commit into from
Sep 20, 2024
Merged

Enable simplex meshes #132

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 21 additions & 13 deletions examples/benchmarks_3D.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,6 @@
#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/mapping_fe.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
Expand Down Expand Up @@ -47,7 +50,7 @@ class AgglomerationBenchmark
make_grid();

Triangulation<dim> tria;
MappingQ<dim> mapping;
MappingFE<dim> mapping;
std::unique_ptr<AgglomerationHandler<dim>> ah;

public:
Expand All @@ -74,7 +77,7 @@ AgglomerationBenchmark<dim>::AgglomerationBenchmark(
const PartitionerType &partitioner_type,
const unsigned int extraction_level,
const unsigned int n_subdomains)
: mapping(1)
: mapping(FE_SimplexDGP<dim>(1))
, grid_type(grid_type)
, partitioner_type(partitioner_type)
, extraction_level(extraction_level)
Expand Down Expand Up @@ -105,9 +108,11 @@ AgglomerationBenchmark<dim>::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
Expand Down Expand Up @@ -158,8 +163,8 @@ AgglomerationBenchmark<dim>::make_grid()


// const auto tree = pack_rtree<bgi::rstar<max_elem_per_node>>(boxes);
auto start = std::chrono::system_clock::now();
const auto tree = pack_rtree<bgi::rstar<max_elem_per_node>>(boxes);
auto start = std::chrono::system_clock::now();
auto tree = pack_rtree<bgi::rstar<max_elem_per_node>>(boxes);
std::cout << "Total number of available levels: " << n_levels(tree)
<< std::endl;

Expand All @@ -169,9 +174,9 @@ AgglomerationBenchmark<dim>::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<vec<>>
CellsAgglomerator<dim, decltype(tree)> 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)
Expand Down Expand Up @@ -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,
Expand All @@ -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();
}
}
232 changes: 232 additions & 0 deletions examples/repairing.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,232 @@
#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/mapping_fe.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>

#include <deal.II/numerics/data_out.h>

#include <agglomeration_handler.h>
#include <poly_utils.h>


// Check that the R-tree based agglomeration works also in 3D.


enum class GridType
{
grid_generator,
unstructured
};

enum class PartitionerType
{
metis,
rtree
};


template <int dim>
class Test
{
private:
void
make_grid();
void
setup_agglomeration();

Triangulation<dim> tria;
MappingFE<dim> mapping;
FE_DGQ<dim> dg_fe;
std::unique_ptr<AgglomerationHandler<dim>> ah;
// no hanging node in DG discretization, we define an AffineConstraints object
// so we can use the distribute_local_to_global() directly.
AffineConstraints<double> constraints;
Vector<double> system_rhs;
std::unique_ptr<GridTools::Cache<dim>> cached_tria;
std::unique_ptr<const Function<dim>> 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 <int dim>
Test<dim>::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<dim>(fe_degree))
, dg_fe(fe_degree)
, grid_type(grid_type)
, partitioner_type(partitioner_type)
, extraction_level(extraction_level)
, n_subdomains(n_subdomains)
{}

template <int dim>
void
Test<dim>::make_grid()
{
GridIn<dim> 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<GridTools::Cache<dim>>(tria, mapping);
}

template <int dim>
void
Test<dim>::setup_agglomeration()

{
ah = std::make_unique<AgglomerationHandler<dim>>(*cached_tria);

std::vector<types::global_cell_index> idxs_to_be_agglomerated = {0, 1, 2, 7};

std::vector<typename Triangulation<2>::active_cell_iterator>
cells_to_be_agglomerated;
PolyUtils::collect_cells_for_agglomeration(tria,
idxs_to_be_agglomerated,
cells_to_be_agglomerated);

std::vector<types::global_cell_index> idxs_to_be_agglomerated2 = {4, 5, 9};

std::vector<typename Triangulation<2>::active_cell_iterator>
cells_to_be_agglomerated2;
PolyUtils::collect_cells_for_agglomeration(tria,
idxs_to_be_agglomerated2,
cells_to_be_agglomerated2);

std::vector<types::global_cell_index> idxs_to_be_agglomerated3 = {3, 6};

std::vector<typename Triangulation<2>::active_cell_iterator>
cells_to_be_agglomerated3;
PolyUtils::collect_cells_for_agglomeration(tria,
idxs_to_be_agglomerated3,
cells_to_be_agglomerated3);

std::vector<types::global_cell_index> idxs_to_be_agglomerated4 = {8,
9,
10,
15};

std::vector<typename Triangulation<2>::active_cell_iterator>
cells_to_be_agglomerated4;
PolyUtils::collect_cells_for_agglomeration(tria,
idxs_to_be_agglomerated4,
cells_to_be_agglomerated4);

std::vector<types::global_cell_index> idxs_to_be_agglomerated5 = {11,
12,
13,
14};

std::vector<typename Triangulation<2>::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<dim> data_out;
data_out.attach_dof_handler(ah->agglo_dh);
// data_out.attach_dof_handler(ah->output_dh);

Vector<float> 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<dim>::type_cell_data);
data_out.build_patches(mapping);
data_out.write_vtu(output);
}
std::cout << "########### Done ###########" << std::endl;
}

template <int dim>
void
Test<dim>::run()
{
make_grid();
setup_agglomeration();
}

int
main()
{
{
Test<2> test{GridType::grid_generator,
PartitionerType::rtree,
4 /*extraction_level*/};
test.run();
}



return 0;
}
20 changes: 17 additions & 3 deletions include/agglomeration_handler.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_fe_field.h>
Expand Down Expand Up @@ -56,9 +57,6 @@
using namespace dealii;

// Forward declarations
template <int dim>
class FinitElement;

template <int dim, int spacedim>
class AgglomerationHandler;

Expand Down Expand Up @@ -298,6 +296,9 @@ class AgglomerationHandler : public Subscriptor
inline const FiniteElement<dim, spacedim> &
get_fe() const;

inline const Mapping<dim> &
get_mapping() const;

inline const std::vector<BoundingBox<dim>> &
get_local_bboxes() const;

Expand Down Expand Up @@ -831,6 +832,10 @@ class AgglomerationHandler : public Subscriptor
// Map the master cell index with the polytope index
std::map<types::global_cell_index, types::global_cell_index> master2polygon;


std::vector<typename Triangulation<dim>::active_cell_iterator>
master_disconnected;

// Dummy FiniteElement objects needed only to generate quadratures

/**
Expand Down Expand Up @@ -883,6 +888,15 @@ AgglomerationHandler<dim, spacedim>::get_fe() const



template <int dim, int spacedim>
inline const Mapping<dim> &
AgglomerationHandler<dim, spacedim>::get_mapping() const
{
return *mapping;
}



template <int dim, int spacedim>
inline const Triangulation<dim, spacedim> &
AgglomerationHandler<dim, spacedim>::get_triangulation() const
Expand Down
Loading
Loading