From c2a3eb17826af821e4a1fa05d13aa6644e138565 Mon Sep 17 00:00:00 2001 From: heister Date: Mon, 5 Mar 2012 16:42:52 +0000 Subject: [PATCH] New: GridGenerator::quarter_hyper_shell() in 3d. git-svn-id: https://svn.dealii.org/trunk@25210 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 4 + deal.II/include/deal.II/grid/grid_generator.h | 61 ++++---- deal.II/source/grid/grid_generator.cc | 136 +++++++++++++++++- 3 files changed, 169 insertions(+), 32 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 177f59bbe2..92fc9e9ffa 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -145,6 +145,10 @@ enabled due to a missing include file in file

Specific improvements

    +
  1. New: GridGenerator::quarter_hyper_shell() in 3d. +
    +(Thomas Geenen, 2012/03/05) +
  2. New: DataOut::write_vtu_in_parallel(). This routine uses MPI I/O to write a big vtu file in parallel.
    diff --git a/deal.II/include/deal.II/grid/grid_generator.h b/deal.II/include/deal.II/grid/grid_generator.h index 6f8d623cce..6406061e86 100644 --- a/deal.II/include/deal.II/grid/grid_generator.h +++ b/deal.II/include/deal.II/grid/grid_generator.h @@ -691,33 +691,28 @@ class GridGenerator /** * Produce a quarter hyper-shell, - * i.e. the space between two - * circles in two space - * dimensions and the region - * between two spheres in 3d, - * with given inner and outer - * radius and a given number of - * elements for this initial - * triangulation. All components are - * restricted to be positive, so the - * opening angle is 90 degrees. - * - * If the number of - * initial cells is zero (as is - * the default), then it is - * computed adaptively such that - * the resulting elements have - * the least aspect ratio. - * - * If colorize is set to true, the - * inner, outer, left, and right - * boundary get indicator 0, 1, 2, - * and 3, respectively. Otherwise - * all indicators are set to 0. + * i.e. the space between two circles in + * two space dimensions and the region + * between two spheres in 3d, with given + * inner and outer radius and a given + * number of elements for this initial + * triangulation. + * + * If the number of initial cells is zero + * (as is the default), then it is + * computed adaptively such that the + * resulting elements have the least + * aspect ratio in 2d. + * + * If colorize is set to true, the inner, + * outer, left, and right boundary get + * indicator 0, 1, 2, and 3 in 2d, + * respectively. Otherwise all indicators + * are set to 0. In 3d indicator 2 is at + * the face x=0, 3 at y=0, 4 at z=0. * * @note The triangulation needs to be - * void upon calling this - * function. Only implemented in 2d so far. + * void upon calling this function. */ template static void quarter_hyper_shell (Triangulation &tria, @@ -988,6 +983,22 @@ class GridGenerator const double inner_radius, const double outer_radius); + + /** + * Assign boundary number zero the inner + * shell boundary, one to the outer shell + * boundary, two to the face with x=0, + * three to the face with y=0, four to + * the face with z=0. + */ + template + static + void + colorize_quarter_hyper_shell(Triangulation & tria, + const Point& center, + const double inner_radius, + const double outer_radius); + /** * Solve the Laplace equation for * @p laplace_transformation diff --git a/deal.II/source/grid/grid_generator.cc b/deal.II/source/grid/grid_generator.cc index 73f6e36614..745ff90461 100644 --- a/deal.II/source/grid/grid_generator.cc +++ b/deal.II/source/grid/grid_generator.cc @@ -2841,16 +2841,138 @@ GridGenerator::half_hyper_shell (Triangulation<3>& tria, } +// Implementation for 3D only +template <> +void +GridGenerator::colorize_quarter_hyper_shell(Triangulation<3> & tria, + const Point<3>& center, + const double inner_radius, + const double outer_radius) +{ + + if (tria.n_cells() != 3) + AssertThrow (false, ExcNotImplemented()); + double middle = (outer_radius-inner_radius)/2e0 + inner_radius; + double eps = 1e-3*middle; + Triangulation<3>::cell_iterator cell = tria.begin(); + + for (;cell!=tria.end();++cell) + for(unsigned int f=0; f::faces_per_cell; ++f) + { + if(!cell->face(f)->at_boundary()) + continue; + + double radius = cell->face(f)->center().norm() - center.norm(); + if (std::fabs(cell->face(f)->center()(0)) < eps ) // x = 0 set boundary 2 + { + cell->face(f)->set_boundary_indicator(2); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) + cell->face(f)->line(j)->set_boundary_indicator(2); + } + else if (std::fabs(cell->face(f)->center()(1)) < eps) // y = 0 set boundary 3 + { + cell->face(f)->set_boundary_indicator(3); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) + cell->face(f)->line(j)->set_boundary_indicator(3); + } + else if (std::fabs(cell->face(f)->center()(2)) < eps ) // z = 0 set boundary 4 + { + cell->face(f)->set_boundary_indicator(4); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) + cell->face(f)->line(j)->set_boundary_indicator(4); + } + else if (radius < middle) // inner radius set boundary 0 + { + cell->face(f)->set_boundary_indicator(0); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps) + cell->face(f)->line(j)->set_boundary_indicator(0); + } + else if (radius > middle) // outer radius set boundary 1 + { + cell->face(f)->set_boundary_indicator(1); + for (unsigned int j=0;j::lines_per_face;++j) + if(cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps) + cell->face(f)->line(j)->set_boundary_indicator(1); + } + else + AssertThrow (false, ExcInternalError()); + } +} + + + +// Implementation for 3D only template <> -void GridGenerator::quarter_hyper_shell (Triangulation<3> &, - const Point<3>&, - const double, - const double, - const unsigned int, - const bool) +void GridGenerator::quarter_hyper_shell (Triangulation<3> & tria, + const Point<3>& center, + const double inner_radius, + const double outer_radius, + const unsigned int n, + const bool colorize) { - Assert (false, ExcNotImplemented()); + Assert ((inner_radius > 0) && (inner_radius < outer_radius), + ExcInvalidRadii ()); + if (n == 0 || n == 3) + { + const double a = inner_radius*std::sqrt(2.0)/2e0; + const double b = outer_radius*std::sqrt(2.0)/2e0; + const double c = a*std::sqrt(3.0)/2e0; + const double d = b*std::sqrt(3.0)/2e0; + const double e = outer_radius/2e0; + const double f = outer_radius*std::sqrt(3.0)/2e0; + const double g = f*std::sqrt(3.0)/2e0; + const double h = inner_radius/2e0; + + std::vector > vertices; + + vertices.push_back (center+Point<3>( 0, inner_radius, 0)); //0 + vertices.push_back (center+Point<3>( a, a, 0)); //1 + vertices.push_back (center+Point<3>( b, b, 0)); //2 + vertices.push_back (center+Point<3>( 0, outer_radius, 0)); //3 + vertices.push_back (center+Point<3>( 0, a , a)); //4 + vertices.push_back (center+Point<3>( c, c , h)); //5 + vertices.push_back (center+Point<3>( d, d , e)); //6 + vertices.push_back (center+Point<3>( 0, b , b)); //7 + vertices.push_back (center+Point<3>( inner_radius, 0 , 0)); //8 + vertices.push_back (center+Point<3>( outer_radius, 0 , 0)); //9 + vertices.push_back (center+Point<3>( a, 0 , a)); //10 + vertices.push_back (center+Point<3>( b, 0 , b)); //11 + vertices.push_back (center+Point<3>( 0, 0 , inner_radius)); //12 + vertices.push_back (center+Point<3>( 0, 0 , outer_radius)); //13 + + const int cell_vertices[3][8] = { + {0, 1, 3, 2, 4, 5, 7, 6}, + {1, 8, 2, 9, 5, 10, 6, 11}, + {4, 5, 7, 6, 12, 10, 13, 11}, + }; + std::vector > cells(n); + + for (unsigned int i=0; i