const Point<3> ¢er = Point<3>(),
const double radius = 1.);
+ /**
+ * Create a cylinder around the
+ * x-axis. The cylinder extends
+ * from @p{x=-half_length} to
+ * @p{x=+half_length} and its
+ * projection into the
+ * @p{yz}-plane is a circle of
+ * radius @p{radius}.
+ *
+ * The boundaries are colored
+ * according to the following
+ * scheme: 0 for the hull of the
+ * cylinder, 1 for the left hand
+ * face and 2 for the right hand
+ * face.
+ */
+ static void cylinder (Triangulation<3> &tria,
+ const double radius = 1.,
+ const double half_length = 1.);
+
+ /**
+ * Projection of the
+ * three-dimensional cylinder
+ * into the @p{xy}-plane.
+ * Therefore, this is simply a square.
+ *
+ * Coloring is like in 3D.
+ */
+ static void cylinder (Triangulation<2> &tria,
+ const double radius = 1.,
+ const double half_length = 1.);
+
+ /**
+ * Non-implemented dummy for compilation purposes.
+ */
+ static void cylinder (Triangulation<1> &tria,
+ const double radius,
+ const double half_length);
+
+
+
/**
* Initialize the given
* triangulation with a hyper-L
#include <grid/tria_boundary.h>
+/**
+ * Boundary object for the hull of a cylinder. In three dimensions,
+ * points are projected on a circular tube along the @p{x}-axis. The
+ * radius of the tube can be set. Similar to @ref{HyperBallBoundary},
+ * new points are projected by dividing the straight line between the
+ * old two points and adjusting the radius in the @p{yz}-plane.
+ *
+ * This class was developed to be used in conjunction with the
+ * @p{cylinder} function of @ref{GridGenerator}. It should be used for
+ * the hull of the cylinder only (boundary indicator 0).
+ *
+ * This class is derived from @ref{StraightBoundary} rather than from
+ * @ref{Boundary}, which would seem natural, since this way we can use the
+ * @ref{StraightBoundary}@p{<dim>::in_between(neighbors)} function.
+ *
+ * @author Guido Kanschat, 2001
+ */
+template <int dim>
+class CylinderBoundary : public StraightBoundary<dim>
+{
+ public:
+ /**
+ * Constructor
+ */
+ CylinderBoundary (const double radius = 1.0);
+
+ /**
+ * Refer to the general documentation of
+ * this class and the documentation of the
+ * base class.
+ */
+ virtual Point<dim>
+ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
+
+ /**
+ * Refer to the general documentation of
+ * this class and the documentation of the
+ * base class.
+ */
+ virtual Point<dim>
+ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
+
+ /**
+ * Refer to the general
+ * documentation of this class
+ * and the documentation of the
+ * base class.
+ *
+ * Calls
+ * @p{get_intermediate_points_between_points}.
+ */
+ virtual void
+ get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
+ typename std::vector<Point<dim> > &points) const;
+
+ /**
+ * Refer to the general
+ * documentation of this class
+ * and the documentation of the
+ * base class.
+ *
+ * Only implemented for @p{dim=3}
+ * and for @p{points.size()==1}.
+ */
+ virtual void
+ get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
+ typename std::vector<Point<dim> > &points) const;
+
+ /**
+ * Compute the normals to the
+ * boundary at the vertices of
+ * the given face.
+ *
+ * Refer to the general
+ * documentation of this class
+ * and the documentation of the
+ * base class.
+ */
+ virtual void
+ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+ typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const;
+
+ /**
+ * Return the radius of the cylinder.
+ */
+ double get_radius () const;
+
+ /**
+ * Exception. Thrown by the
+ * @p{get_radius} if the
+ * @p{compute_radius_automatically},
+ * see below, flag is set true.
+ */
+ DeclException0 (ExcRadiusNotSet);
+
+
+ protected:
+ /**
+ * Radius of the cylinder.
+ */
+ const double radius;
+
+ private:
+
+ /**
+ * Called by
+ * @p{get_intermediate_points_on_line}
+ * and by
+ * @p{get_intermediate_points_on_quad}.
+ *
+ * Refer to the general
+ * documentation of
+ * @p{get_intermediate_points_on_line}
+ * in the documentation of the
+ * base class.
+ */
+ void get_intermediate_points_between_points (const Point<dim> &p0, const Point<dim> &p1,
+ typename std::vector<Point<dim> > &points) const;
+};
+
+
+
/**
* Specialisation of @ref{Boundary}<dim>, which places the new point on
* the boundary of a ball in arbitrary dimension. It works by projecting
cell->face(3)->set_boundary_indicator (1);
cell->face(4)->set_boundary_indicator (5);
cell->face(5)->set_boundary_indicator (0);
+ break;
default:
Assert(false, ExcNotImplemented());
};
+void
+GridGenerator::cylinder (Triangulation<2> &,
+ const double,
+ const double)
+{
+ Assert (false, ExcInternalError());
+}
+
+
+
+void GridGenerator::cylinder (Triangulation<1> &,
+ const double,
+ const double)
+{
+ Assert (false, ExcInternalError());
+};
+
+
+
void GridGenerator::hyper_shell (Triangulation<1> &,
const Point<1> &,
const double,
+void
+GridGenerator::cylinder (Triangulation<2> &tria,
+ const double radius,
+ const double half_length)
+{
+ Point<2> p1 (-half_length, -radius);
+ Point<2> p2 (half_length, radius);
+
+ hyper_rectangle(tria, p1, p2, true);
+
+ Triangulation<2>::face_iterator f = tria.begin_face();
+ Triangulation<2>::face_iterator end = tria.end_face();
+ while (f != end)
+ {
+ switch (f->boundary_indicator())
+ {
+ case 0:
+ f->set_boundary_indicator(1);
+ break;
+ case 1:
+ f->set_boundary_indicator(2);
+ break;
+ default:
+ f->set_boundary_indicator(0);
+ break;
+ }
+ ++f;
+ }
+}
+
+
+
void
GridGenerator::half_hyper_shell (Triangulation<2> &tria,
const Point<2> ¢er,
+void
+GridGenerator::cylinder (Triangulation<3> &tria,
+ const double radius,
+ const double half_length)
+{
+ // Copy the base from hyper_ball<2>
+ // and transform it to yz
+ const double d = radius/std::sqrt(2.0);
+ const double a = d/(1+std::sqrt(2.0));
+ const Point<3> vertices[24] = {
+ Point<3>(-half_length, -d,-d),
+ Point<3>(-half_length, d,-d),
+ Point<3>(-half_length, -a,-a),
+ Point<3>(-half_length, a,-a),
+ Point<3>(-half_length, -a, a),
+ Point<3>(-half_length, a, a),
+ Point<3>(-half_length, -d, d),
+ Point<3>(-half_length, d, d),
+ Point<3>(0, -d,-d),
+ Point<3>(0, d,-d),
+ Point<3>(0, -a,-a),
+ Point<3>(0, a,-a),
+ Point<3>(0, -a, a),
+ Point<3>(0, a, a),
+ Point<3>(0, -d, d),
+ Point<3>(0, d, d),
+ Point<3>(half_length, -d,-d),
+ Point<3>(half_length, d,-d),
+ Point<3>(half_length, -a,-a),
+ Point<3>(half_length, a,-a),
+ Point<3>(half_length, -a, a),
+ Point<3>(half_length, a, a),
+ Point<3>(half_length, -d, d),
+ Point<3>(half_length, d, d),
+ };
+
+ int cell_vertices[10][8] = {
+ {0,1,3,2,8,9,11,10},
+ {0,2,4,6,8,10,12,14},
+ {2,3,5,4,10,11,13,12},
+ {1,7,5,3,9,15,13,11},
+ {6,4,5,7,14,12,13,15}
+ };
+ 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<CellData<3> > 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<Point<3> >(&vertices[0], &vertices[24]),
+ cells,
+ SubCellData()); // no boundary information
+
+ Triangulation<3>::cell_iterator cell = tria.begin();
+ Triangulation<3>::cell_iterator end = tria.end();
+
+ while (cell != end)
+ {
+ if (cell->face(0)->boundary_indicator() != 255)
+ cell->face(0)->set_boundary_indicator(1);
+ if (cell->face(1)->boundary_indicator() != 255)
+ cell->face(1)->set_boundary_indicator(2);
+ ++cell;
+ }
+};
+
+
+
void GridGenerator::hyper_shell (Triangulation<3> &,
const Point<3> &,
const double ,
void GridOut::write_dx (const Triangulation<dim> &tria,
std::ostream &out)
{
- AssertThrow (out, ExcIO());
-
- // get the positions of the
- // vertices and whether they are
- // used.
- const std::vector<Point<dim> > &vertices = tria.get_vertices();
- const std::vector<bool> &vertex_used = tria.get_used_vertices();
-
- const unsigned int n_vertices = tria.n_used_vertices();
-
- typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
- const typename Triangulation<dim>::active_cell_iterator endc=tria.end();
-
- // start with ucd data
- out << n_vertices << ' '
- << tria.n_active_cells() + (ucd_flags.write_faces ?
- n_boundary_faces(tria) :
- 0)
- << " 0 0 0" // no data
- << std::endl;
-
- // actually write the vertices.
- // note that we shall number them
- // with first index 1 instead of 0
- for (unsigned int i=0; i<vertices.size(); ++i)
- if (vertex_used[i])
- {
- out << i+1 // vertex index
- << " "
- << vertices[i];
- for (unsigned int d=dim+1; d<=3; ++d)
- out << " 0"; // fill with zeroes
- out << std::endl;
- };
-
- // write cells. Enumerate cells
- // consecutively, starting with 1
- unsigned int cell_index=1;
- for (cell=tria.begin_active();
- cell!=endc; ++cell, ++cell_index)
- {
- out << cell_index << ' '
- << static_cast<unsigned int>(cell->material_id())
- << " ";
- switch (dim)
- {
- case 1: out << "line "; break;
- case 2: out << "quad "; break;
- case 3: out << "hex "; break;
- default:
- Assert (false, ExcNotImplemented());
- };
-
- // it follows a list of the
- // vertices of each cell. in 1d
- // this is simply a list of the
- // two vertices, in 2d its counter
- // clockwise, as usual in this
- // library. in 3d, the same applies
- // (special thanks to AVS for
- // numbering their vertices in a
- // way compatible to deal.II!)
- //
- // technical reference:
- // AVS Developer's Guide, Release 4,
- // May, 1992, p. E6
- //
- // note: vertex numbers are 1-base
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- out << cell->vertex_index(vertex)+1 << ' ';
- out << std::endl;
- };
-
- // write faces with non-zero boundary
- // indicator
- if (ucd_flags.write_faces)
- write_ucd_faces (tria, cell_index, out);
-
+ Assert(false, ExcNotImplemented());
AssertThrow (out, ExcIO());
};
#include <cmath>
+template <int dim>
+CylinderBoundary<dim>::CylinderBoundary (const double radius) :
+ radius(radius)
+{};
+
+
+
+template <int dim>
+Point<dim>
+CylinderBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+{
+ Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
+
+ // project to boundary
+ if (dim>=3)
+ middle *= radius / std::sqrt(middle.square()-middle(0)*middle(0));
+ return middle;
+};
+
+
+template <int dim>
+Point<dim>
+CylinderBoundary<dim>::
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
+{
+ Point<dim> middle = StraightBoundary<dim>::get_new_point_on_quad (quad);
+
+ // project to boundary
+ if (dim>=3)
+ middle *= radius / std::sqrt(middle.square()-middle(0)*middle(0));
+ return middle;
+};
+
+
+template <int dim>
+void
+CylinderBoundary<dim>::get_intermediate_points_on_line (
+ const typename Triangulation<dim>::line_iterator &line,
+ typename std::vector<Point<dim> > &points) const
+{
+ if (points.size()==1)
+ points[0]=get_new_point_on_line(line);
+ else
+ get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
+}
+
+
+template <int dim>
+void
+CylinderBoundary<dim>::get_intermediate_points_between_points (
+ const Point<dim> &v0, const Point<dim> &v1,
+ typename std::vector<Point<dim> > &points) const
+{
+ const unsigned int n=points.size();
+ Assert(n>0, ExcInternalError());
+ // Do a simple linear interpolation
+ // followed by projection. If this
+ // is not sufficient, try to
+ // understand the sophisticated
+ // code in HyperBall later.
+ Point<dim> ds = v1-v0;
+ ds /= n+1;
+
+ for (unsigned int i=0; i<n; ++i)
+ {
+ if (i==0)
+ points[i] = v0+ds;
+ else
+ points[i] = points[i-1]+ds;
+
+ points[i] *= radius / std::sqrt(points[i].square()-points[i](0)*points[i](0));
+ }
+}
+
+
+#if deal_II_dimension == 3
+
+template <>
+void
+CylinderBoundary<3>::get_intermediate_points_on_quad (
+ const Triangulation<3>::quad_iterator &quad,
+ std::vector<Point<3> > &points) const
+{
+ unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
+ Assert(points.size()==m*m, ExcInternalError());
+
+ std::vector<Point<3> > lp3(m);
+ std::vector<Point<3> > lp1(m);
+
+ get_intermediate_points_on_line(quad->line(3), lp3);
+ get_intermediate_points_on_line(quad->line(1), lp1);
+
+ std::vector<Point<3> > lps(m);
+ for (unsigned int i=0; i<m; ++i)
+ {
+ get_intermediate_points_between_points(lp3[i], lp1[i], lps);
+
+ for (unsigned int j=0; j<m; ++j)
+ points[i*m+j]=lps[j];
+ }
+}
+
+#endif
+
+
+template <int dim>
+void
+CylinderBoundary<dim>::get_intermediate_points_on_quad (
+ const typename Triangulation<dim>::quad_iterator &,
+ typename std::vector<Point<dim> > &) const
+{
+ Assert(false, Boundary<dim>::ExcFunctionNotUseful(dim));
+}
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+void
+CylinderBoundary<1>::
+get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+ Boundary<1>::FaceVertexNormals &) const
+{
+ Assert (false, Boundary<1>::ExcFunctionNotUseful(1));
+};
+
+#endif
+
+
+template <int dim>
+void
+CylinderBoundary<dim>::
+get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+ typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
+{
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+ {
+ face_vertex_normals[vertex] = face->vertex(vertex);
+ face_vertex_normals[vertex][0] = 0.;
+ }
+};
+
+
+
+template <int dim>
+double
+CylinderBoundary<dim>::get_radius () const
+{
+ return radius;
+};
+
+
+//======================================================================//
+
template <int dim>
HyperBallBoundary<dim>::HyperBallBoundary (const Point<dim> p,
const double radius) :
// explicit instantiations
template class HyperBallBoundary<deal_II_dimension>;
+template class CylinderBoundary<deal_II_dimension>;
template class HalfHyperBallBoundary<deal_II_dimension>;
template class HyperShellBoundary<deal_II_dimension>;
template class HalfHyperShellBoundary<deal_II_dimension>;