From 5fac5406aaf634579a835e23553d3b0dd6e94fde Mon Sep 17 00:00:00 2001 From: heltai Date: Thu, 22 Aug 2013 03:05:52 +0000 Subject: [PATCH] Default manifold is set to boundary_id when that one is set and manifold_id is invalid. Added implementation of CylindricalManifold. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@30403 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-1/step-1.cc | 11 +- deal.II/include/deal.II/grid/manifold_lib.h | 289 ++++++++---------- .../deal.II/grid/tria_accessor.templates.h | 12 +- deal.II/source/grid/manifold_lib.cc | 77 +++++ deal.II/source/grid/manifold_lib.inst.in | 16 +- 5 files changed, 213 insertions(+), 192 deletions(-) diff --git a/deal.II/examples/step-1/step-1.cc b/deal.II/examples/step-1/step-1.cc index 42def746f6..724dd45e95 100644 --- a/deal.II/examples/step-1/step-1.cc +++ b/deal.II/examples/step-1/step-1.cc @@ -105,16 +105,7 @@ void second_grid () Triangulation<2>::active_cell_iterator cell = triangulation.begin_active(), - endc = triangulation.end(); - for(; cell != endc; ++cell) - cell->set_all_manifold_ids(0); - - - // for(unsigned int f=0; f::faces_per_cell; ++f) - // if(cell->face(f)->at_boundary()) - // cell->face(f)->set_manifold_id(0); - - + endc = triangulation.end(); // By default, the triangulation assumes that all boundaries are straight // and given by the cells of the coarse grid (which we just created). It diff --git a/deal.II/include/deal.II/grid/manifold_lib.h b/deal.II/include/deal.II/grid/manifold_lib.h index 97d203e8fc..e426dfa67f 100644 --- a/deal.II/include/deal.II/grid/manifold_lib.h +++ b/deal.II/include/deal.II/grid/manifold_lib.h @@ -25,190 +25,143 @@ DEAL_II_NAMESPACE_OPEN /** - * Manifold description for a polar space coordinate system. Given any - * two points in space, tranform them in polar coordinate centered at - * center, compute intermediate points in polar coordinate systems, - * and returns its space representation. + * Manifold description for a polar space coordinate system. Given a + * set of points in space, this class computes their weighted average + * and the weighted average of their distance from the center. The + * average of the points is then shifted in the radial direction, to + * match the weighted average of the distances. + * + * You can use this Manifold object to describe any sphere, circle, + * hypersphere or hyperdisc in two or three dimensions, both as a + * co-dimension one manifold descriptor or as co-dimension zero + * manifold descriptor. + * + * @ingroup manifold + * + * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007, Luca Heltai, + * 2013 */ template class PolarManifold : public FlatManifold { public: + /** + * The Constructor takes the + * center of the polar + * coordinates system. + */ PolarManifold(const Point center = Point()); + /** + * Given a set of points in + * space, this function computes + * their weighted average and the + * weighted average of their + * distance from the center. The + * average of the points is then + * shifted in the radial + * direction, to match the + * weighted average of the + * distances. + */ virtual Point get_new_point (const std::vector > &surrounding_points, const std::vector &weights) const; private: - Point center; + /** + * The center of the polar + * coordinate system. + */ + const Point center; +}; + +/** + * Manifold description for a cylindrical space coordinate system. In + * three dimensions, points are projected on a cylindrical tube along + * the x-, y- or z-axis (when using the + * first constructor of this class), or an arbitrarily oriented + * cylinder described by the direction of its axis and a point located + * on the axis. Similar to PolarManifold, new points are projected by + * taking the weighetd averages of the old points and adjusting the + * radius from the axis. + * + * This class was developed to be used in conjunction with the @p + * cylinder function of GridGenerator. Its use is discussed in detail + * in the results section of step-49. + * + * @ingroup manifold + * + * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007, Luca Heltai, + * 2013 + */ + +template +class CylinderManifold : public FlatManifold +{ + public: + /** + * Constructor. Using default + * values for the constructor + * arguments yields a cylindrical + * manifold tube along the x-axis + * (axis=0). Choose + * axis=1 or + * axis=2 for a + * cylindrical manifold along the + * y- or z-axis, respectively. + */ + CylinderManifold(const int axis=0); + + /** + * Constructor. If constructed + * with this constructor, the + * boundary described is a + * cylinder with an axis that + * points in direction #direction + * and goes through the given + * #point_on_axis. The direction + * may be arbitrarily scaled, and + * the given point may be any + * point on the axis. + */ + CylinderManifold(const Point &direction, + const Point &point_on_axis); + + /** + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual Point + get_new_point (const std::vector > &surrounding_points, + const std::vector &weights) const; + + + private: + /** + * Direction of the cylindrical + * coordinate system. + */ + const Point direction; + + /** + * One point on the axis of the + * cylindrical coordinate systems + */ + const Point point_on_axis; + + /** + * Returns the axis vector + * associated with the given axis + * number. + */ + static Point get_axis_vector (const unsigned int axis); }; -// -// /** -// * Boundary object for the hull of a cylinder. In three dimensions, -// * points are projected on a circular tube along the x-, -// * y- or z-axis (when using the first constructor of -// * this class), or an arbitrarily oriented cylinder described by the -// * direction of its axis and a point located on the axis. The radius -// * of the tube can be given independently. Similar to -// * HyperBallBoundary, new points are projected by dividing the -// * straight line between the old two points and adjusting the radius -// * from the axis. -// * -// * This class was developed to be used in conjunction with the -// * @p cylinder function of GridGenerator. It should be used for -// * the hull of the cylinder only (boundary indicator 0). Its use is -// * discussed in detail in the results section of step-49. -// * -// * This class is derived from StraightBoundary rather than from -// * Boundary, which would seem natural, since this way we can use the -// * StraightBoundary::in_between() function. -// * -// * @ingroup boundary -// * -// * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007 -// */ -// template -// class CylinderBoundary : public StraightBoundary -// { -// public: -// /** -// * Constructor. Using default values for the constructor arguments yields a -// * circular tube along the x-axis (axis=0). Choose axis=1 -// * or axis=2 for a tube along the y- or z-axis, respectively. -// */ -// CylinderBoundary (const double radius = 1.0, -// const unsigned int axis = 0); -// -// /** -// * Constructor. If constructed -// * with this constructor, the -// * boundary described is a -// * cylinder with an axis that -// * points in direction #direction -// * and goes through the given -// * #point_on_axis. The direction -// * may be arbitrarily scaled, and -// * the given point may be any -// * point on the axis. -// */ -// CylinderBoundary (const double radius, -// const Point &direction, -// const Point &point_on_axis); -// -// /** -// * Refer to the general documentation of -// * this class and the documentation of the -// * base class. -// */ -// virtual Point -// get_new_point_on_line (const typename Triangulation::line_iterator &line) const; -// -// /** -// * Refer to the general documentation of -// * this class and the documentation of the -// * base class. -// */ -// virtual Point -// get_new_point_on_quad (const typename Triangulation::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::line_iterator &line, -// std::vector > &points) const; -// -// /** -// * Refer to the general -// * documentation of this class -// * and the documentation of the -// * base class. -// * -// * Only implemented for dim=3 -// * and for points.size()==1. -// */ -// virtual void -// get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, -// std::vector > &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::face_iterator &face, -// typename Boundary::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; -// -// /** -// * The direction vector of the axis. -// */ -// const Point direction; -// -// /** -// * An arbitrary point on the axis. -// */ -// const Point point_on_axis; -// -// 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 &p0, const Point &p1, -// std::vector > &points) const; -// -// /** -// * Given a number for the axis, -// * return a vector that denotes -// * this direction. -// */ -// static Point get_axis_vector (const unsigned int axis); -// }; // // // diff --git a/deal.II/include/deal.II/grid/tria_accessor.templates.h b/deal.II/include/deal.II/grid/tria_accessor.templates.h index e7d0d57325..14fade0c7b 100644 --- a/deal.II/include/deal.II/grid/tria_accessor.templates.h +++ b/deal.II/include/deal.II/grid/tria_accessor.templates.h @@ -1939,8 +1939,6 @@ template const Boundary & TriaAccessor::get_boundary () const { - // This assert should be removed in - // future releases... Assert (structdimused(), TriaAccessorExceptions::ExcCellNotUsed()); // Get the default (manifold_id) @@ -1968,7 +1966,17 @@ TriaAccessor::get_manifold () const Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed()); // Get the default (manifold_id) types::manifold_id mi = this->objects().manifold_id[this->present_index]; + + // In case this is not valid, check + // the boundary id, after having + // casted it to a manifold id + if(mi == numbers::invalid_manifold_id) + mi = static_cast + (structdim < dim ? + this->objects().boundary_or_material_id[this->present_index].boundary_id: + this->objects().boundary_or_material_id[this->present_index].material_id); + // Now recover the manifold itself return this->tria->get_manifold(mi); } diff --git a/deal.II/source/grid/manifold_lib.cc b/deal.II/source/grid/manifold_lib.cc index 5df357c341..a7b5ba4a32 100644 --- a/deal.II/source/grid/manifold_lib.cc +++ b/deal.II/source/grid/manifold_lib.cc @@ -49,8 +49,85 @@ get_new_point (const std::vector > &surrounding_points, return p; } +// ============================================================ +// CylinderManifold +// ============================================================ +template +CylinderManifold::CylinderManifold(const int axis): + direction (get_axis_vector (axis)), + point_on_axis (Point()) +{} + +template +CylinderManifold::CylinderManifold(const Point &direction, + const Point &point_on_axis): + direction (direction), + point_on_axis (point_on_axis) +{} + + + +template +Point +CylinderManifold::get_axis_vector (const unsigned int axis) +{ + + Assert (axis < spacedim, ExcIndexRange (axis, 0, spacedim)); + + Point axis_vector; + + axis_vector[axis] = 1; + + return axis_vector; +} + + +template +Point +CylinderManifold:: +get_new_point (const std::vector > &surrounding_points, + const std::vector &weights) const +{ + Assert(surrounding_points.size() == weights.size(), + ExcDimensionMismatch(surrounding_points.size(), weights.size())); + + // compute a proposed new point + Point middle = FlatManifold::get_new_point(surrounding_points, weights); + + double radius = 0; + Point on_plane; + + for(unsigned int i=0; i vector_from_axis = (middle-point_on_axis) - + ((middle-point_on_axis) * direction) * direction; + + // scale it to the desired length + // and put everything back + // together, unless we have a point + // on the axis + if (vector_from_axis.norm() <= 1e-10 * middle.norm()) + return middle; + + else + return (vector_from_axis / vector_from_axis.norm() * radius + + ((middle-point_on_axis) * direction) * direction + + point_on_axis); +} + // explicit instantiations #include "manifold_lib.inst" diff --git a/deal.II/source/grid/manifold_lib.inst.in b/deal.II/source/grid/manifold_lib.inst.in index 2a5ddfc978..1191c0d1fe 100644 --- a/deal.II/source/grid/manifold_lib.inst.in +++ b/deal.II/source/grid/manifold_lib.inst.in @@ -17,19 +17,11 @@ for (deal_II_dimension : DIMENSIONS) { -// template class PolarManifold; +#if deal_II_dimension != 1 + template class PolarManifold; + template class CylinderManifold; +#endif // template class ConeBoundary; -// template class HyperBallBoundary; -// template class HalfHyperBallBoundary; -// template class HyperShellBoundary; -// template class HalfHyperShellBoundary; -// -// #if deal_II_dimension != 3 -// template class HyperBallBoundary; -// #endif -// #if deal_II_dimension == 3 -// template class CylinderBoundary; -// #endif } -- 2.39.5