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

Add new Mapping based on agglomerates #133

Merged
merged 4 commits into from
Sep 29, 2024
Merged

Add new Mapping based on agglomerates #133

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
13b0f1a
Add new class MappingBox, implementing bbox mapping on general (hex o…
fdrmrc Sep 12, 2024
c1e7de4
Replace euler mapping with box mapping
fdrmrc Sep 16, 2024
36b3e38
Remove leftovers
fdrmrc Sep 16, 2024
dbff07d
Improve allocation of qpoints inside reinit()
fdrmrc Sep 17, 2024
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
682 changes: 135 additions & 547 deletions examples/minimal_SIP.cc
View file
Open in desktop

Large diffs are not rendered by default.

226 changes: 71 additions & 155 deletions examples/poisson.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,30 +12,27 @@


#include <deal.II/fe/fe_dgp.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/grid/reference_cell.h>

#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_matrix.h>

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

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

#include <algorithm>
#include <chrono>


#define HEX TRUE

struct ConvergenceInfo
{
Expand Down Expand Up @@ -73,17 +70,6 @@ struct ConvergenceInfo



template <typename T>
constexpr T
constexpr_pow(T num, unsigned int pow)
{
return (pow >= sizeof(unsigned int) * 8) ? 0 :
pow == 0 ? 1 :
num * constexpr_pow(num, pow - 1);
}



enum class GridType
{
grid_generator, // hyper_cube or hyper_ball
Expand Down Expand Up @@ -421,9 +407,14 @@ class Poisson
output_results();


Triangulation<dim> tria;
MappingQ<dim> mapping;
FE_DGQ<dim> dg_fe;
Triangulation<dim> tria;
#ifdef HEX
MappingQ1<dim> mapping;
FE_DGQ<dim> dg_fe;
#else
MappingFE<dim> mapping;
FE_SimplexDGP<dim> dg_fe;
#endif
std::unique_ptr<AgglomerationHandler<dim>> ah;
AffineConstraints<double> constraints;
SparsityPattern sparsity;
Expand Down Expand Up @@ -470,7 +461,12 @@ Poisson<dim>::Poisson(const GridType &grid_type,
const unsigned int extraction_level,
const unsigned int n_subdomains,
const unsigned int fe_degree)
: mapping(1)
:
#ifdef HEX
mapping()
#else
mapping(FE_SimplexP<dim>{1})
#endif
, dg_fe(fe_degree)
, grid_type(grid_type)
, partitioner_type(partitioner_type)
Expand Down Expand Up @@ -502,30 +498,42 @@ Poisson<dim>::make_grid()
if constexpr (dim == 2)
{
grid_in.attach_triangulation(tria);
#ifdef HEX
std::ifstream gmsh_file(
"../../meshes/t3.msh"); // unstructured square made by triangles
#else
std::ifstream gmsh_file(
"../../meshes/t3.msh"); // unstructured square [0,1]^2
"../../meshes/square_simplex_coarser.msh"); // unstructured square
// made by triangles
#endif
grid_in.read_msh(gmsh_file);
tria.refine_global(5); // 4
tria.refine_global(2); // 4
}
else if constexpr (dim == 3)
{
// {
// grid_in.attach_triangulation(tria);
// std::ifstream gmsh_file("../../meshes/piston_3.inp"); //
// piston
// mesh grid_in.read_abaqus(gmsh_file); tria.refine_global(1);
// }
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);
#ifdef HEX
std::ifstream filename("../../meshes/piston_3.inp"); // piston mesh
grid_in.read_abaqus(filename);
#else
std::ifstream filename(
"../../meshes/gray_level_image1.vtk"); // liver or brain domain
grid_in.read_vtk(filename);

#endif
}
}
else
{
#ifdef HEX
GridGenerator::hyper_cube(tria, 0., 1.);
tria.refine_global(8);
tria.refine_global(4);
#else
Triangulation<dim> tria_hex;
GridGenerator::hyper_cube(tria_hex, 0., 1.);
tria_hex.refine_global(4);
GridGenerator::convert_hypercube_to_simplex_mesh(tria_hex, tria);
#endif
}
std::cout << "Size of tria: " << tria.n_active_cells() << std::endl;
cached_tria = std::make_unique<GridTools::Cache<dim>>(tria, mapping);
Expand All @@ -550,7 +558,7 @@ Poisson<dim>::make_grid()

namespace bgi = boost::geometry::index;
static constexpr unsigned int max_elem_per_node =
constexpr_pow(2, dim); // 2^dim
PolyUtils::constexpr_pow(2, dim); // 2^dim
std::vector<std::pair<BoundingBox<dim>,
typename Triangulation<dim>::active_cell_iterator>>
boxes(tria.n_active_cells());
Expand Down Expand Up @@ -693,42 +701,6 @@ Poisson<dim>::setup_agglomeration()
data_out.build_patches(mapping);
data_out.write_vtu(output);
}



/*
#ifdef AGGLO_DEBUG
{
for (const auto &cell : ah->agglomeration_cell_iterators())
{
std::cout << "Cell with idx: " << cell->active_cell_index()
<< std::endl;
unsigned int n_agglomerated_faces_per_cell = ah->n_faces(cell);
std::cout << "Number of faces for the agglomeration: "
<< n_agglomerated_faces_per_cell << std::endl;
for (unsigned int f = 0; f < n_agglomerated_faces_per_cell; ++f)
{
std::cout << "Agglomerated face with idx: " << f << std::endl;
const auto &[local_face_idx, neigh, local_face_idx_out, dummy] =
ah->get_agglomerated_connectivity()[{cell, f}];
{
std::cout << "Face idx: " << local_face_idx << std::endl;
if (neigh.state() == IteratorState::valid)
{
std::cout << "Neighbor idx: " <<
neigh->active_cell_index()
<< std::endl;
}

std::cout << "Face idx from outside: " << local_face_idx_out
<< std::endl;
}
std::cout << std::endl;
}
}
}
#endif
*/
}


Expand All @@ -741,15 +713,24 @@ Poisson<dim>::assemble_system()
solution.reinit(ah->n_dofs());
system_rhs.reinit(ah->n_dofs());

const unsigned int quadrature_degree = 2 * dg_fe.get_degree() + 1;
const unsigned int face_quadrature_degree = 2 * dg_fe.get_degree() + 1;
const unsigned int quadrature_degree = dg_fe.get_degree() + 1;
const unsigned int face_quadrature_degree = dg_fe.get_degree() + 1;
#ifdef HEX
ah->initialize_fe_values(QGauss<dim>(quadrature_degree),
update_gradients | update_JxW_values |
update_quadrature_points | update_JxW_values |
update_values,
QGauss<dim - 1>(face_quadrature_degree));
#else
ah->initialize_fe_values(QGaussSimplex<dim>(quadrature_degree),
update_gradients | update_JxW_values |
update_quadrature_points | update_JxW_values |
update_values,
QGaussSimplex<dim - 1>(face_quadrature_degree));
#endif

const unsigned int dofs_per_cell = ah->n_dofs_per_cell();
std::cout << "DoFs per cell: " << dofs_per_cell << std::endl;

FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
Expand Down Expand Up @@ -911,7 +892,7 @@ Poisson<dim>::assemble_system()
++i)
{
double d = (points0[i] - points1[i]).norm();
Assert(
AssertThrow(
d < 1e-15,
ExcMessage(
"Face qpoints at the interface do not match!"));
Expand All @@ -930,11 +911,6 @@ Poisson<dim>::assemble_system()

const auto &normals = fe_faces0.get_normal_vectors();

// const double penalty =
// penalty_constant /
// PolyUtils::compute_h_orthogonal(f,
// polygon_boundary_vertices,
// normals[0]);
const double penalty =
penalty_constant / std::fabs(polytope->diameter());

Expand Down Expand Up @@ -1027,12 +1003,6 @@ Poisson<dim>::assemble_system()



void
output_double_number(double input, const std::string &text)
{
std::cout << text << input << std::endl;
}

template <int dim>
void
Poisson<dim>::solve()
Expand All @@ -1048,27 +1018,6 @@ template <int dim>
void
Poisson<dim>::output_results()
{
{
const std::string filename = "agglomerated_Poisson.vtu";
std::ofstream output(filename);

DataOut<dim> data_out;
data_out.attach_dof_handler(ah->agglo_dh);
data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
data_out.build_patches(/**(ah->euler_mapping)*/ mapping);
data_out.write_vtu(output);
}
{
const std::string filename = "agglomerated_Poisson_test.vtu";
std::ofstream output(filename);

DataOut<dim> data_out;
data_out.attach_dof_handler(ah->agglo_dh);
data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
data_out.build_patches(*(ah->euler_mapping), 3);
data_out.write_vtu(output);
}

{
std::string partitioner;
if (partitioner_type == PartitionerType::metis)
Expand All @@ -1078,13 +1027,16 @@ Poisson<dim>::output_results()
else
partitioner = "no_partitioning";

const std::string filename = "interpolated_solution_" + partitioner + "_" +
const std::string filename = "interpolated_solution" + partitioner + "_" +
std::to_string(n_subdomains) + ".vtu";
std::ofstream output(filename);

DataOut<dim> data_out;
Vector<double> interpolated_solution;
PolyUtils::interpolate_to_fine_grid(*ah, interpolated_solution, solution);
PolyUtils::interpolate_to_fine_grid(*ah,
interpolated_solution,
solution,
true /*no_the_fly*/);
data_out.attach_dof_handler(ah->output_dh);
data_out.add_data_vector(interpolated_solution,
"u",
Expand Down Expand Up @@ -1120,53 +1072,6 @@ Poisson<dim>::output_results()
semih1_err = errors[1];
std::cout << "Error (L2): " << l2_err << std::endl;
std::cout << "Error (H1): " << semih1_err << std::endl;

// Old version
#ifdef FALSE
{
// L2 error
Vector<float> difference_per_cell(tria.n_active_cells());
VectorTools::integrate_difference(mapping,
ah->output_dh,
interpolated_solution,
*analytical_solution,
difference_per_cell,
QGauss<dim>(dg_fe.degree),
VectorTools::L2_norm);

// Plot the error per cell

data_out.add_data_vector(difference_per_cell,
"error_per_cell",
DataOut<dim>::type_cell_data);

const double L2_error =
VectorTools::compute_global_error(tria,
difference_per_cell,
VectorTools::L2_norm);

std::cout << "L2 error:" << L2_error << std::endl;
}

{
// H1 seminorm
Vector<float> difference_per_cell(tria.n_active_cells());
VectorTools::integrate_difference(mapping,
ah->output_dh,
interpolated_solution,
*analytical_solution,
difference_per_cell,
QGauss<dim>(dg_fe.degree + 1),
VectorTools::H1_seminorm);

const double H1_seminorm =
VectorTools::compute_global_error(tria,
difference_per_cell,
VectorTools::H1_seminorm);

std::cout << "H1 seminorm:" << H1_seminorm << std::endl;
}
#endif
}
}

Expand Down Expand Up @@ -1196,7 +1101,14 @@ Poisson<dim>::run()
{
make_grid();
setup_agglomeration();
auto start = std::chrono::high_resolution_clock::now();
assemble_system();
auto stop = std::chrono::high_resolution_clock::now();
auto duration =
std::chrono::duration_cast<std::chrono::microseconds>(stop - start);

std::cout << "Time taken by assemble_system(): " << duration.count() / 1e6
<< " seconds" << std::endl;
solve();
output_results();
}
Expand All @@ -1210,7 +1122,11 @@ main()
ConvergenceInfo convergence_info;
std::cout << "Testing p-convergence" << std::endl;
{
for (unsigned int fe_degree : {1, 2, 3, 4, 5, 6})
#ifdef HEX
for (unsigned int fe_degree : {1, 2, 3, 4})
#else
for (unsigned int fe_degree : {1, 2, 3})
#endif
{
std::cout << "Fe degree: " << fe_degree << std::endl;
Poisson<2> poisson_problem{GridType::unstructured,
Expand Down
Loading
Loading