From f25bbcc4b721d9043349a8111dce28b0ff061ecf Mon Sep 17 00:00:00 2001 From: Bruno Date: Thu, 1 Apr 2021 16:40:22 -0400 Subject: [PATCH] Regular cylinder uses subdivided_cylinder --- source/grid/grid_generator.cc | 117 +--------------------------------- 1 file changed, 1 insertion(+), 116 deletions(-) diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 3f2fb6382e..160681e6ef 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -5461,122 +5461,7 @@ namespace GridGenerator const double radius, const double half_length) { - // Copy the base from hyper_ball<3> - // and transform it to yz - const double d = radius / std::sqrt(2.0); - const double a = d / (1 + std::sqrt(2.0)); - Point<3> vertices[24] = { - Point<3>(-d, -half_length, -d), - Point<3>(d, -half_length, -d), - Point<3>(-a, -half_length, -a), - Point<3>(a, -half_length, -a), - Point<3>(-a, -half_length, a), - Point<3>(a, -half_length, a), - Point<3>(-d, -half_length, d), - Point<3>(d, -half_length, d), - Point<3>(-d, 0, -d), - Point<3>(d, 0, -d), - Point<3>(-a, 0, -a), - Point<3>(a, 0, -a), - Point<3>(-a, 0, a), - Point<3>(a, 0, a), - Point<3>(-d, 0, d), - Point<3>(d, 0, d), - Point<3>(-d, half_length, -d), - Point<3>(d, half_length, -d), - Point<3>(-a, half_length, -a), - Point<3>(a, half_length, -a), - Point<3>(-a, half_length, a), - Point<3>(a, half_length, a), - Point<3>(-d, half_length, d), - Point<3>(d, half_length, d), - }; - // Turn cylinder such that y->x - for (auto &vertex : vertices) - { - const double h = vertex(1); - vertex(1) = -vertex(0); - vertex(0) = h; - } - - int cell_vertices[10][8] = {{0, 1, 8, 9, 2, 3, 10, 11}, - {0, 2, 8, 10, 6, 4, 14, 12}, - {2, 3, 10, 11, 4, 5, 12, 13}, - {1, 7, 9, 15, 3, 5, 11, 13}, - {6, 4, 14, 12, 7, 5, 15, 13}}; - for (unsigned int i = 0; i < 5; ++i) - for (unsigned int j = 0; j < 8; ++j) - cell_vertices[i + 5][j] = cell_vertices[i][j] + 8; - - std::vector> cells(10, CellData<3>()); - - for (unsigned int i = 0; i < 10; ++i) - { - for (unsigned int j = 0; j < 8; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; - } - - tria.create_triangulation(std::vector>(std::begin(vertices), - std::end(vertices)), - cells, - SubCellData()); // no boundary information - - // set boundary indicators for the - // faces at the ends to 1 and 2, - // respectively. note that we also - // have to deal with those lines - // that are purely in the interior - // of the ends. we determine whether - // an edge is purely in the - // interior if one of its vertices - // is at coordinates '+-a' as set - // above - Triangulation<3>::cell_iterator cell = tria.begin(); - Triangulation<3>::cell_iterator end = tria.end(); - - tria.set_all_manifold_ids_on_boundary(0); - - for (; cell != end; ++cell) - for (unsigned int i : GeometryInfo<3>::face_indices()) - if (cell->at_boundary(i)) - { - if (cell->face(i)->center()(0) > half_length - 1.e-5) - { - cell->face(i)->set_boundary_id(2); - cell->face(i)->set_manifold_id(numbers::flat_manifold_id); - - for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face; - ++e) - if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a)) - { - cell->face(i)->line(e)->set_boundary_id(2); - cell->face(i)->line(e)->set_manifold_id( - numbers::flat_manifold_id); - } - } - else if (cell->face(i)->center()(0) < -half_length + 1.e-5) - { - cell->face(i)->set_boundary_id(1); - cell->face(i)->set_manifold_id(numbers::flat_manifold_id); - - for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face; - ++e) - if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a)) - { - cell->face(i)->line(e)->set_boundary_id(1); - cell->face(i)->line(e)->set_manifold_id( - numbers::flat_manifold_id); - } - } - } - tria.set_manifold(0, CylindricalManifold<3>()); + subdivided_cylinder(tria, 2, radius, half_length); } // Implementation for 3D only -- 2.39.5