From 8d3afa82e175a535ef8bb5ea465ad4b78e3871b1 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 9 Sep 2011 02:23:06 +0000 Subject: [PATCH] Patch by Daniel Castanon Quiroz: Implement generating a torus and describing the boundary. git-svn-id: https://svn.dealii.org/trunk@24288 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 6 +- deal.II/include/deal.II/grid/grid_generator.h | 49 ++- .../include/deal.II/grid/tria_boundary_lib.h | 116 ++++++- deal.II/source/grid/grid_generator.cc | 155 ++++++++- deal.II/source/grid/tria_boundary_lib.cc | 298 ++++++++++++++++++ 5 files changed, 593 insertions(+), 31 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 6adb15686d..2055fe1710 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -283,13 +283,17 @@ and DoF handlers embedded in higher dimensional space have been added.

Specific improvements

    +
  1. New: The GridGenerator::torus function and TorusBoundary class can +create and describe the surface of a torus. +
    +(Daniel Castanon Quiroz, 2011/09/08) +
  2. Fixed: FEFieldFunction class was not thread-safe because it keeps a cache on the side that was invalidated when different threads kept pouncing on it. This is now fixed.
    (Patrick Sodré, 2011/09/07) -
  3. New: There is now a function GridTools::volume computing the volume of a triangulation.
    diff --git a/deal.II/include/deal.II/grid/grid_generator.h b/deal.II/include/deal.II/grid/grid_generator.h index 0e5bc83357..9981b9563a 100644 --- a/deal.II/include/deal.II/grid/grid_generator.h +++ b/deal.II/include/deal.II/grid/grid_generator.h @@ -47,9 +47,6 @@ template class SparseMatrix; * * * @ingroup grid - * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat, Stefan - * Nauber, Joerg Weimar, Yaqi Wang, Luca Heltai, 1998, 1999, 2000, 2001, 2002, - * 2003, 2006, 2007, 2008, 2009, 2010, 2011. */ class GridGenerator { @@ -673,11 +670,11 @@ class GridGenerator * 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 + * and 3, respectively. Otherwise * all indicators are set to 0. * * @note The triangulation needs to be @@ -692,7 +689,7 @@ class GridGenerator const unsigned int n_cells = 0, const bool colorize = false); - + /** * Produce a quarter hyper-shell, * i.e. the space between two @@ -730,7 +727,7 @@ class GridGenerator const double outer_radius, const unsigned int n_cells = 0, const bool colorize = false); - + /** * Produce a domain that is the space * between two cylinders in 3d, with @@ -761,6 +758,30 @@ class GridGenerator const unsigned int n_radial_cells = 0, const unsigned int n_axial_cells = 0); + + + /** + * Produce the surface meshing of the + * torus. The axis of the torus is the + * $y$-axis while the plane of the torus + * is the $x$-$z$ plane. The boundary of + * this object can be described by the + * TorusBoundary class. + * + * @param R The radius of the circle, + * which forms the middle line of the + * torus containing the loop of + * cells. Must be greater than @p r. + * + * @param r The inner radius of the + * torus. + */ + + static void torus (Triangulation<2,3>& tria, + const double R, + const double r); + + /** * This class produces a square * on the xy-plane with a @@ -826,8 +847,8 @@ class GridGenerator * the third parameter. Previous * content of @p result will be * deleted. - * - * This function is most often used + * + * This function is most often used * to compose meshes for more * complicated geometries if the * geometry can be composed of @@ -835,7 +856,7 @@ class GridGenerator * exist to generate coarse meshes. * For example, the channel mesh used * in step-35 could in principle be - * created using a mesh created by the + * created using a mesh created by the * GridGenerator::hyper_cube_with_cylindrical_hole * function and several rectangles, * and merging them using the current @@ -847,11 +868,11 @@ class GridGenerator * individual mesh building blocks are * GridTools::transform, GridTools::rotate, * and GridTools::scale). - * + * * @note The two input triangulations * must be coarse meshes that have - * no refined cells. - * + * no refined cells. + * * @note The function copies the material ids * of the cells of the two input * triangulations into the output @@ -863,7 +884,7 @@ class GridGenerator * then you will currently have to set * boundary indicators again by hand * in the output triangulation. - * + * * @note For a related operation * on refined meshes when both * meshes are derived from the diff --git a/deal.II/include/deal.II/grid/tria_boundary_lib.h b/deal.II/include/deal.II/grid/tria_boundary_lib.h index eb5f68f804..42ef1c8681 100644 --- a/deal.II/include/deal.II/grid/tria_boundary_lib.h +++ b/deal.II/include/deal.II/grid/tria_boundary_lib.h @@ -388,7 +388,7 @@ class HyperBallBoundary : public StraightBoundary * this class and the documentation of the * base class. */ - virtual + virtual Point get_new_point_on_line (const typename Triangulation::line_iterator &line) const; @@ -397,7 +397,7 @@ class HyperBallBoundary : public StraightBoundary * this class and the documentation of the * base class. */ - virtual + virtual Point get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; @@ -410,7 +410,7 @@ class HyperBallBoundary : public StraightBoundary * Calls * @p get_intermediate_points_between_points. */ - virtual + virtual void get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, std::vector > &points) const; @@ -424,7 +424,7 @@ class HyperBallBoundary : public StraightBoundary * Only implemented for dim=3 * and for points.size()==1. */ - virtual + virtual void get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, std::vector > &points) const; @@ -453,7 +453,7 @@ class HyperBallBoundary : public StraightBoundary * and the documentation of the * base class. */ - virtual + virtual void get_normals_at_vertices (const typename Triangulation::face_iterator &face, typename Boundary::FaceVertexNormals &face_vertex_normals) const; @@ -467,7 +467,7 @@ class HyperBallBoundary : public StraightBoundary /** * Return the radius of the ball. */ - double + double get_radius () const; /** @@ -734,6 +734,109 @@ class HalfHyperShellBoundary : public HyperShellBoundary }; +/** + * Class describing the boundary of the torus. The axis of the torus is the + * $y$-axis while the plane of the torus is the $x$-$z$ plane. A torus of this + * kind can be generated by GridGenerator::torus. + * + * This class is only implemented for + * dim=2,spacedim=3, that is, just the surface. + */ +template +class TorusBoundary : public Boundary +{ + public: + /** + * Constructor.R has to be + * greater than r. + */ + TorusBoundary (const double R, const double r); + +//Boundary Refinenment Functions + /** + * Construct a new point on a line. + */ + virtual Point + get_new_point_on_line (const typename Triangulation::line_iterator &line) const; + + /** + * Construct a new point on a quad. + */ + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; + + /** + * Construct a new points on a line. + */ + virtual void get_intermediate_points_on_line ( + const typename Triangulation< dim, spacedim >::line_iterator & line, + std::vector< Point< spacedim > > & points) const ; + + /** + * Construct a new points on a quad. + */ + virtual void get_intermediate_points_on_quad ( + const typename Triangulation< dim, spacedim >::quad_iterator & quad, + std::vector< Point< spacedim > > & points ) const ; + + /** + * Get the normal from cartesian + * coordinates. This normal does not have + * unit length. + */ + virtual void get_normals_at_vertices ( + const typename Triangulation< dim, spacedim >::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const ; + + private: + //Handy functions + /** + * Function that corrects the value and + * sign of angle, that is, given + * angle=tan(abs(y/x)); if + * (y > 0) && (x < 0) then + * correct_angle = Pi - angle, + * etc. + */ + + double get_correct_angle(const double angle,const double x,const double y) const; + + /** + * Get the cartesian coordinates of the + * Torus, i.e., from (theta,phi) + * to (x,y,z). + */ + Point get_real_coord(const Point& surfP) const ; + + /** + * Get the surface coordinates of the + * Torus, i.e., from (x,y,z) to + * (theta,phi). + */ + Point get_surf_coord(const Point& p) const ; + + /** + * Get the normal from surface + * coordinates. This normal does not have + * unit length. + */ + Point get_surf_norm_from_sp(const Point& surfP) const ; + + /** + * Get the normal from cartesian + * coordinates. This normal does not have + * unit length. + */ + Point get_surf_norm(const Point& p) const ; + + /** + * Inner and outer radii of the shell. + */ + const double R; + const double r; +}; + + /* -------------- declaration of explicit specializations ------------- */ @@ -787,6 +890,7 @@ HalfHyperShellBoundary<1>:: get_normals_at_vertices (const Triangulation<1,1>::face_iterator &, Boundary<1,1>::FaceVertexNormals &) const; + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/grid/grid_generator.cc b/deal.II/source/grid/grid_generator.cc index 515f92b1f2..1e4bd92f11 100644 --- a/deal.II/source/grid/grid_generator.cc +++ b/deal.II/source/grid/grid_generator.cc @@ -232,6 +232,141 @@ GridGenerator::moebius ( +void +GridGenerator::torus (Triangulation<2,3>& tria, + const double R, + const double r) +{ + Assert (R>r, ExcMessage("Outer radius must be greater than inner radius.")); + + const unsigned int dim=2; + const unsigned int spacedim=3; + std::vector > vertices (16); + + vertices[0]=Point(R-r,0,0); + vertices[1]=Point(R,-r,0); + vertices[2]=Point(R+r,0,0); + vertices[3]=Point(R, r,0); + vertices[4]=Point(0,0,R-r); + vertices[5]=Point(0,-r,R); + vertices[6]=Point(0,0,R+r); + vertices[7]=Point(0,r,R); + vertices[8]=Point(-(R-r),0,0); + vertices[9]=Point(-R,-r,0); + vertices[10]=Point(-(R+r),0,0); + vertices[11]=Point(-R, r,0); + vertices[12]=Point(0,0,-(R-r)); + vertices[13]=Point(0,-r,-R); + vertices[14]=Point(0,0,-(R+r)); + vertices[15]=Point(0,r,-R); + + std::vector > cells (16); + //Right Hand Orientation + cells[0].vertices[0] = 0; + cells[0].vertices[1] = 4; + cells[0].vertices[2] = 7; + cells[0].vertices[3] = 3; + cells[0].material_id = 0; + + cells[1].vertices[0] = 1; + cells[1].vertices[1] = 5; + cells[1].vertices[2] = 4; + cells[1].vertices[3] = 0; + cells[1].material_id = 0; + + cells[2].vertices[0] = 2; + cells[2].vertices[1] = 6; + cells[2].vertices[2] = 5; + cells[2].vertices[3] = 1; + cells[2].material_id = 0; + + cells[3].vertices[0] = 3; + cells[3].vertices[1] = 7; + cells[3].vertices[2] = 6; + cells[3].vertices[3] = 2; + cells[3].material_id = 0; + + cells[4].vertices[0] = 4; + cells[4].vertices[1] = 8; + cells[4].vertices[2] = 11; + cells[4].vertices[3] = 7; + cells[4].material_id = 0; + + cells[5].vertices[0] = 5; + cells[5].vertices[1] = 9; + cells[5].vertices[2] = 8; + cells[5].vertices[3] = 4; + cells[5].material_id = 0; + + cells[6].vertices[0] = 6; + cells[6].vertices[1] = 10; + cells[6].vertices[2] = 9; + cells[6].vertices[3] = 5; + cells[6].material_id = 0; + + cells[7].vertices[0] = 7; + cells[7].vertices[1] = 11; + cells[7].vertices[2] = 10; + cells[7].vertices[3] = 6; + cells[7].material_id = 0; + + cells[8].vertices[0] = 8; + cells[8].vertices[1] = 12; + cells[8].vertices[2] = 15; + cells[8].vertices[3] = 11; + cells[8].material_id = 0; + + cells[9].vertices[0] = 9; + cells[9].vertices[1] = 13; + cells[9].vertices[2] = 12; + cells[9].vertices[3] = 8; + cells[9].material_id = 0; + + cells[10].vertices[0] = 10; + cells[10].vertices[1] = 14; + cells[10].vertices[2] = 13; + cells[10].vertices[3] = 9; + cells[10].material_id = 0; + + cells[11].vertices[0] = 11; + cells[11].vertices[1] = 15; + cells[11].vertices[2] = 14; + cells[11].vertices[3] = 10; + cells[11].material_id = 0; + + cells[12].vertices[0] = 12; + cells[12].vertices[1] = 0; + cells[12].vertices[2] = 3; + cells[12].vertices[3] = 15; + cells[12].material_id = 0; + + cells[13].vertices[0] = 13; + cells[13].vertices[1] = 1; + cells[13].vertices[2] = 0; + cells[13].vertices[3] = 12; + cells[13].material_id = 0; + + cells[14].vertices[0] = 14; + cells[14].vertices[1] = 2; + cells[14].vertices[2] = 1; + cells[14].vertices[3] = 13; + cells[14].material_id = 0; + + cells[15].vertices[0] = 15; + cells[15].vertices[1] = 3; + cells[15].vertices[2] = 2; + cells[15].vertices[3] = 14; + cells[15].material_id = 0; + + // Must call this to be able to create a + // correct triangulation in dealii, read + // GridReordering<> doc + GridReordering::reorder_cells (cells); + tria.create_triangulation_compatibility (vertices, cells, SubCellData()); +} + + + // Implementation for 2D only template<> @@ -1697,7 +1832,7 @@ GridGenerator::half_hyper_shell (Triangulation<2> &tria, }; tria.create_triangulation (vertices, cells, SubCellData()); - + if (colorize) { Triangulation<2>::cell_iterator cell = tria.begin(); @@ -1706,7 +1841,7 @@ GridGenerator::half_hyper_shell (Triangulation<2> &tria, cell->face(2)->set_boundary_indicator(1); } tria.begin()->face(0)->set_boundary_indicator(3); - + tria.last()->face(1)->set_boundary_indicator(2); } } @@ -1776,8 +1911,8 @@ void GridGenerator::quarter_hyper_shell (Triangulation<2> &tria, cells[i].material_id = 0; }; - tria.create_triangulation (vertices, cells, SubCellData()); - + tria.create_triangulation (vertices, cells, SubCellData()); + if (colorize) { Triangulation<2>::cell_iterator cell = tria.begin(); @@ -1786,7 +1921,7 @@ void GridGenerator::quarter_hyper_shell (Triangulation<2> &tria, cell->face(2)->set_boundary_indicator(1); } tria.begin()->face(0)->set_boundary_indicator(3); - + tria.last()->face(1)->set_boundary_indicator(2); } } @@ -2813,13 +2948,13 @@ merge_triangulations (const Triangulation &triangulation_1, ExcMessage ("The input triangulations must be coarse meshes.")); Assert (triangulation_2.n_levels() == 1, ExcMessage ("The input triangulations must be coarse meshes.")); - + // get the union of the set of vertices std::vector > vertices = triangulation_1.get_vertices(); vertices.insert (vertices.end(), triangulation_2.get_vertices().begin(), triangulation_2.get_vertices().end()); - + // now form the union of the set of cells std::vector > cells; cells.reserve (triangulation_1.n_cells() + triangulation_2.n_cells()); @@ -2832,8 +2967,8 @@ merge_triangulations (const Triangulation &triangulation_1, this_cell.material_id = cell->material_id(); cells.push_back (this_cell); } - - // now do the same for the other other mesh. note that we have to + + // now do the same for the other other mesh. note that we have to // translate the vertex indices for (typename Triangulation::cell_iterator cell = triangulation_2.begin(); cell != triangulation_2.end(); ++cell) @@ -2844,7 +2979,7 @@ merge_triangulations (const Triangulation &triangulation_1, this_cell.material_id = cell->material_id(); cells.push_back (this_cell); } - + // throw out duplicated vertices from the two meshes // and create the triangulation SubCellData subcell_data; diff --git a/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/source/grid/tria_boundary_lib.cc index 2cf45eac64..88074b4c19 100644 --- a/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/source/grid/tria_boundary_lib.cc @@ -1220,6 +1220,304 @@ get_normals_at_vertices (const typename Triangulation::face_iterator &face, + +template +TorusBoundary::TorusBoundary (const double R__, + const double r__) + : + R(R__), + r(r__) +{ + Assert (false, ExcNotImplemented()); +} + + + +template <> +TorusBoundary<2,3>::TorusBoundary (const double R__, + const double r__) + : + R(R__), + r(r__) +{ + Assert (R>r, ExcMessage("Outer radius must be greater than inner radius.")); +} + + + +template +double +TorusBoundary::get_correct_angle(const double angle, + const double x, + const double y) const +{ + if(y>=0) + { + if (x >=0) + return angle; + + return numbers::PI-angle; + } + + if (x <=0) + return numbers::PI+angle; + + return 2.0*numbers::PI-angle; +} + + + +template <> +Point<3> +TorusBoundary<2,3>::get_real_coord (const Point<2>& surfP) const +{ + const double theta=surfP(0); + const double phi=surfP(1); + + return Point<3> ((R+r*std::cos(phi))*std::cos(theta), + r*std::sin(phi), + (R+r*std::cos(phi))*std::sin(theta)); +} + + + +template <> +Point<2> +TorusBoundary<2,3>::get_surf_coord(const Point<3>& p) const +{ + const double phi=std::asin(std::abs(p(1))/r); + const double Rr_2=p(0)*p(0)+p(2)*p(2); + + Point<2> surfP; + surfP(1)=get_correct_angle(phi,Rr_2-R*R,p(1));//phi + + if(std::abs(p(0))<1.E-5) + { + if(p(2)>=0) + surfP(0) = numbers::PI*0.5; + else + surfP(0) = -numbers::PI*0.5; + } + else + { + const double theta = std::atan(std::abs(p(2)/p(0))); + surfP(0)= get_correct_angle(theta,p(0),p(2)); + } + + return surfP; +} + + + +template <> +Point<3> +TorusBoundary<2,3>::get_new_point_on_line (const typename Triangulation<2,3>::line_iterator &line) const +{ + //Just get the average + Point<2> p0=get_surf_coord(line->vertex(0)); + Point<2> p1=get_surf_coord(line->vertex(1)); + + Point<2> middle(0,0); + + //Take care for periodic conditions, + //For instance phi0= 0, phi1= 3/2*Pi middle has to be 7/4*Pi not 3/4*Pi + //This also works for -Pi/2 + Pi, middle is 5/4*Pi + for(unsigned i=0; i<2;i++) + if(std::abs(p0(i)-p1(i))> numbers::PI) + middle(i)=2*numbers::PI; + + middle+= p0 + p1; + middle*=0.5; + + Point<3> midReal=get_real_coord(middle); + return midReal; +} + + + +template <> +Point<3> +TorusBoundary<2,3>::get_new_point_on_quad (const typename Triangulation<2,3>::quad_iterator &quad) const +{ + //Just get the average + Point<2> p[4]; + + for(unsigned i=0;i<4;i++) + p[i]=get_surf_coord(quad->vertex(i)); + + Point<2> middle(0,0); + + //Take care for periodic conditions, see get_new_point_on_line() above + //For instance phi0= 0, phi1= 3/2*Pi middle has to be 7/4*Pi not 3/4*Pi + //This also works for -Pi/2 + Pi + Pi- Pi/2, middle is 5/4*Pi + for(unsigned i=0;i<2;i++) + for(unsigned j=1;j<4;j++){ + if(std::abs(p[0](i)-p[j](i))> numbers::PI){ + middle(i)+=2*numbers::PI; + } + } + + for(unsigned i=0;i<4;i++) + middle+=p[i]; + + middle*= 0.25; + + return get_real_coord(middle); +} + + + +//Normal field without unit length +template <> +Point<3> +TorusBoundary<2,3>:: get_surf_norm_from_sp(const Point<2>& surfP) const +{ + + Point<3> n; + double theta=surfP[0]; + double phi=surfP[1]; + + double f=R+r*std::cos(phi); + + n[0]=r*std::cos(phi)*std::cos(theta)*f; + n[1]=r*std::sin(phi)*f; + n[2]=r*std::sin(theta)*std::cos(phi)*f; + + return n; +} + + + +//Normal field without unit length +template <> +Point<3> +TorusBoundary<2,3>::get_surf_norm(const Point<3>& p) const{ + + Point<2> surfP=get_surf_coord(p); + return get_surf_norm_from_sp(surfP); + +} + + + +template<> +void +TorusBoundary<2,3>:: +get_intermediate_points_on_line (const typename Triangulation<2, 3>::line_iterator & line, + std::vector< Point< 3 > > & points) const +{ + //Almost the same implementation as StraightBoundary<2,3> + unsigned npoints=points.size(); + if(npoints==0) return; + + Point<2> p[2]; + + for(unsigned i=0;i<2;i++) + p[i]=get_surf_coord(line->vertex(i)); + + unsigned offset[2]; + offset[0]=0; + offset[1]=0; + + //Take care for periodic conditions & negative angles, + //see get_new_point_on_line() above + //Because we dont have a symmetric interpolation (just the middle) we need to + //add 2*Pi to each almost zero and negative angles. + for(unsigned i=0;i<2;i++) + for(unsigned j=1;j<2;j++){ + if(std::abs(p[0](i)-p[j](i))> numbers::PI){ + offset[i]++; + break; + } + } + + for(unsigned i=0;i<2;i++) + for(unsigned j=0;j<2;j++) + if(p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles + p[j](i)+=2*numbers::PI*offset[i]; + + + double dx=1.0/(npoints+1); + double x=dx; + Point<2> target; + for(unsigned i=0; i +void +TorusBoundary<2,3>:: +get_intermediate_points_on_quad (const typename Triangulation< 2, 3 >::quad_iterator &quad, + std::vector< Point< 3 > > & points )const +{ + //Almost the same implementation as StraightBoundary<2,3> + const unsigned int n=points.size(), + m=static_cast(std::sqrt(static_cast(n))); + // is n a square number + Assert(m*m==n, ExcInternalError()); + + const double ds=1./(m+1); + double y=ds; + + Point<2> p[4]; + + for(unsigned i=0;i<4;i++) + p[i]=get_surf_coord(quad->vertex(i)); + + Point<2> target; + unsigned offset[2]; + offset[0]=0; + offset[1]=0; + + //Take care for periodic conditions & negative angles, + //see get_new_point_on_line() above + //Because we dont have a symmetric interpolation (just the middle) we need to + //add 2*Pi to each almost zero and negative angles. + for(unsigned i=0;i<2;i++) + for(unsigned j=1;j<4;j++){ + if(std::abs(p[0](i)-p[j](i))> numbers::PI){ + offset[i]++; + break; + } + } + + for(unsigned i=0;i<2;i++) + for(unsigned j=0;j<4;j++) + if(p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles + p[j](i)+=2*numbers::PI*offset[i]; + + for (unsigned i=0; i +void +TorusBoundary<2,3>:: +get_normals_at_vertices (const typename Triangulation<2,3 >::face_iterator &face, + FaceVertexNormals &face_vertex_normals) const +{ + for(unsigned i=0; i::vertices_per_face; i++) + face_vertex_normals[i]=get_surf_norm(face->vertex(i)); +} + + + // explicit instantiations #include "tria_boundary_lib.inst" -- 2.39.5