From 54e649a8d9072694569f337aaa4d3ef4cab36820 Mon Sep 17 00:00:00 2001 From: heltai Date: Tue, 20 Aug 2013 15:35:27 +0000 Subject: [PATCH] Added Manifold interface. The branch compiles, no tests have been performed. Documentation needs to be fixed. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@30358 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/grid/manifold.h | 172 ++++ deal.II/include/deal.II/grid/manifold_lib.h | 900 ++++++++++++++++++ deal.II/include/deal.II/grid/tria.h | 35 +- deal.II/include/deal.II/grid/tria_accessor.h | 13 + .../deal.II/grid/tria_accessor.templates.h | 14 + deal.II/source/grid/CMakeLists.txt | 3 + deal.II/source/grid/manifold.cc | 72 ++ deal.II/source/grid/manifold.inst.in | 26 + deal.II/source/grid/tria.cc | 416 +++++--- deal.II/source/grid/tria.inst.in | 2 +- 10 files changed, 1506 insertions(+), 147 deletions(-) create mode 100644 deal.II/include/deal.II/grid/manifold.h create mode 100644 deal.II/include/deal.II/grid/manifold_lib.h create mode 100644 deal.II/source/grid/manifold.cc create mode 100644 deal.II/source/grid/manifold.inst.in diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h new file mode 100644 index 0000000000..8354314640 --- /dev/null +++ b/deal.II/include/deal.II/grid/manifold.h @@ -0,0 +1,172 @@ +// --------------------------------------------------------------------- +// $Id: tria_boundary.h 30130 2013-07-23 13:01:18Z heltai $ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef __deal2__manifold_h +#define __deal2__manifold_h + + +/*---------------------------- boundary-function.h ---------------------------*/ + +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +template class Triangulation; + + + +/** + * This class is used to represent a geometry to a triangulation. + * When a triangulation creates a new vertex in the + * domain, it determines the new vertex' coordinates through the + * following code (here in two dimensions): + * @code + * ... + * Point<2> new_vertex = manifold.get_new_point(surrounding_points, weights); + * ... + * @endcode + * @p surrounding_points is a vector containing all points which surround + * the vertex to be created, while weights contains a list of the weights to use + * when computing the new location. + * + * For example, for a flat manifold, the new_vertex is computed as: + * @code + * ... + * Point new_vertex; + * + * for(unsigned int i=0; i +class Manifold : public Subscriptor +{ +public: + + /** + * Destructor. Does nothing here, but + * needs to be declared to make it + * virtual. + */ + virtual ~Manifold (); + + /** + * Return the point which shall + * become the new vertex surrounded + * by the given points. + */ + virtual + Point + get_new_point(const std::vector > &surrounding_points, + const std::vector &weights) const = 0; + + + /** + * Given a candidate point and a + * line segment characterized by + * the iterator, return a point + * that lies on the surface + * described by this object. This + * function is used in some mesh + * smoothing algorithms that try + * to move around points in order + * to improve the mesh quality + * but need to ensure that points + * that were on the boundary + * remain on the boundary. + * + * If spacedim==1, then the line + * represented by the line + * iterator is the entire space + * (i.e. it is a cell, not a part + * of the boundary), and the + * returned point equals the + * given input point. + * + * Derived classes do not need to + * implement this function unless + * mesh smoothing algorithms are + * used with a particular + * boundary object. The default + * implementation of this + * function throws an exception + * of type ExcPureFunctionCalled. + */ + virtual + Point + project_to_manifold (const Point &candidate) const; +}; + + + +/** + * Specialization of Manifold, which places the new point + * right into the middle of the given points. The middle is defined + * as the weighted mean of the points. + * + * @ingroup manifold + * + * @author Luca Heltai, 2013 + */ +template +class FlatManifold: public Manifold +{ +public: + /** + * Default constructor. Some + * compilers require this for + * some reasons. + */ + FlatManifold (); + + /** + * Let the new point be the + * arithmetic mean of the two + * vertices of the line. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class for more + * information. + */ + virtual Point + get_new_point(const std::vector > &surrounding_points, + const std::vector &weights) const; +}; + + + +/* -------------- declaration of explicit specializations ------------- */ + +#ifndef DOXYGEN + +#endif // DOXYGEN + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/include/deal.II/grid/manifold_lib.h b/deal.II/include/deal.II/grid/manifold_lib.h new file mode 100644 index 0000000000..7bd1cf5e29 --- /dev/null +++ b/deal.II/include/deal.II/grid/manifold_lib.h @@ -0,0 +1,900 @@ +// --------------------------------------------------------------------- +// $Id: tria_boundary_lib.h 30130 2013-07-23 13:01:18Z heltai $ +// +// Copyright (C) 1999 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef __deal2__manifold_lib_h +#define __deal2__manifold_lib_h + + +#include +#include + +DEAL_II_NAMESPACE_OPEN +// +// /** +// * 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); +// }; +// +// +// +// /** +// * Boundary object for the hull of a (truncated) cone with two +// * different radii at the two ends. If one radius is chosen to be 0 +// * the object describes the boundary of a cone. In three dimensions, +// * points are projected on an arbitrarily oriented (truncated) cone +// * described by the two endpoints and the corresponding radii. 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 is derived from StraightBoundary rather than from +// * Boundary, which would seem natural, since this way we can use the +// * StraightBoundary::in_between() function. +// * +// * As an example of use, consider the following code snippet: +// * @code +// * Triangulation triangulation; +// * GridGenerator::truncated_cone (triangulation); +// * Point p1, p2; +// * p1[0] = -1; +// * p2[0] = 1; +// * const ConeBoundary boundary (1, 0.5, p1, p2); +// * triangulation.set_boundary (0, boundary); +// * triangulation.refine_global (2); +// * @endcode +// * This will produce the following meshes after the two +// * refinements we perform, in 2d and 3d, respectively: +// * +// * @image html cone_2d.png +// * @image html cone_3d.png +// * +// * @author Markus Bürg, 2009 +// */ +// template +// class ConeBoundary : public StraightBoundary +// { +// public: +// /** +// * Constructor. Here the boundary +// * object is constructed. The +// * points x_0 and +// * x_1 describe the +// * starting and ending points of +// * the axis of the (truncated) +// * cone. radius_0 +// * denotes the radius +// * corresponding to x_0 +// * and radius_1 the one +// * corresponding to x_1. +// */ +// ConeBoundary (const double radius_0, +// const double radius_1, +// const Point x_0, +// const Point x_1); +// +// /** +// * Return the radius of the +// * (truncated) cone at given point +// * x on the axis. +// */ +// double get_radius (const Point x) const; +// +// /** +// * 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; +// +// protected: +// /** +// * First radius of the (truncated) +// * cone. +// */ +// const double radius_0; +// +// /** +// * Second radius of the (truncated) +// * cone. +// */ +// const double radius_1; +// +// /** +// * Starting point of the axis. +// */ +// const Point x_0; +// +// /** +// * Ending point of the axis. +// */ +// const Point x_1; +// +// 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; +// }; +// +// +// +// /** +// * Specialization of Boundary, which places the new point on +// * the boundary of a ball in arbitrary dimension. It works by projecting +// * the point in the middle of the old points onto the ball. The middle is +// * defined as the arithmetic mean of the points. +// * +// * The center of the ball and its radius may be given upon construction of +// * an object of this type. They default to the origin and a radius of 1.0. +// * +// * 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 Wolfgang Bangerth, 1998, Ralf Hartmann, 2001 +// */ +// template +// class HyperBallBoundary : public StraightBoundary +// { +// public: +// /** +// * Constructor +// */ +// HyperBallBoundary (const Point p = Point(), +// const double radius = 1.0); +// +// /** +// * 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; +// +// /** +// * Implementation of the function +// * declared in the base class. +// * +// * Refer to the general +// * documentation of this class +// * and the documentation of the +// * base class. +// */ +// virtual +// Tensor<1,spacedim> +// normal_vector (const typename Triangulation::face_iterator &face, +// const Point &p) 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 center of the ball. +// */ +// Point +// get_center () const; +// +// /** +// * Return the radius of the ball. +// */ +// 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: +// +// /** +// * Center point of the hyperball. +// */ +// const Point center; +// +// /** +// * Radius of the hyperball. +// */ +// const double radius; +// +// /** +// * This flag is @p false for +// * this class and for all derived +// * classes that set the radius by +// * the constructor. For example +// * this flag is @p false for the +// * HalfHyperBallBoundary +// * class but it is @p true for +// * the HyperShellBoundary +// * class, for example. The +// * latter class doesn't get its +// * radii by the constructor but +// * need to compute the radii +// * automatically each time one of +// * the virtual functions is +// * called. +// */ +// bool compute_radius_automatically; +// +// 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; +// }; +// +// +// +// /** +// * Variant of HyperBallBoundary which denotes a half hyper ball +// * where the first coordinate is restricted to the range $x>=0$ (or +// * $x>=center(0)$). In two dimensions, this equals the right half +// * circle, in three space dimensions it is a half ball. This class +// * might be useful for computations with rotational symmetry, where +// * one dimension is the radius from the axis of rotation. +// * +// * @ingroup boundary +// * +// * @author Wolfgang Bangerth, 1999, 2001 +// */ +// template +// class HalfHyperBallBoundary : public HyperBallBoundary +// { +// public: +// /** +// * Constructor +// */ +// HalfHyperBallBoundary (const Point p = Point(), +// const double radius = 1.0); +// +// /** +// * Check if on the line x==0, +// * otherwise pass to the base +// * class. +// */ +// virtual Point +// get_new_point_on_line (const typename Triangulation::line_iterator &line) const; +// +// /** +// * Check if on the line x==0, +// * otherwise pass to 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; +// }; +// +// +// +// /** +// * Class describing the boundaries of a hyper shell. Only the center +// * of the two spheres needs to be given, the radii of inner and outer +// * sphere are computed automatically upon calling one of the virtual +// * functions. +// * +// * @ingroup boundary +// * +// * @author Wolfgang Bangerth, 1999 +// */ +// template +// class HyperShellBoundary : public HyperBallBoundary +// { +// public: +// /** +// * Constructor. The center of the +// * spheres defaults to the +// * origin. +// * +// * Calls the constructor of its +// * base @p HyperBallBoundary +// * class with a dummy radius as +// * argument. This radius will be +// * ignored +// */ +// HyperShellBoundary (const Point ¢er = Point()); +// }; +// +// +// +// /** +// * Variant of HyperShellBoundary which denotes a half hyper shell +// * where the first coordinate is restricted to the range $x>=0$ (or +// * $x>=center(0)$). In two dimensions, this equals the right half arc, +// * in three space dimensions it is a half shell. This class might be +// * useful for computations with rotational symmetry, where one +// * dimension is the radius from the axis of rotation. +// * +// * @ingroup boundary +// * +// * @author Wolfgang Bangerth, 2000, 2009 +// */ +// template +// class HalfHyperShellBoundary : public HyperShellBoundary +// { +// public: +// /** +// * Constructor. The center of the +// * spheres defaults to the +// * origin. +// * +// * If the radii are not specified, the +// * class tries to infer them from the +// * location of points on the +// * boundary. This works in 2d, but not in +// * 3d. As a consequence, in 3d these +// * radii must be given. +// */ +// HalfHyperShellBoundary (const Point ¢er = Point(), +// const double inner_radius = -1, +// const double outer_radius = -1); +// +// /** +// * 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; +// +// /** +// * 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; +// +// private: +// /** +// * Inner and outer radii of the shell. +// */ +// const double inner_radius; +// const double outer_radius; +// }; +// +// +// /** +// * 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 ------------- */ +// +// #ifndef DOXYGEN +// +// template <> +// Point<1> +// HyperBallBoundary<1>:: +// get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const; +// template <> +// void +// HyperBallBoundary<1>::get_intermediate_points_on_line ( +// const Triangulation<1>::line_iterator &, +// std::vector > &) const; +// template <> +// void +// HyperBallBoundary<3>::get_intermediate_points_on_quad ( +// const Triangulation<3>::quad_iterator &quad, +// std::vector > &points) const; +// template <> +// void +// HyperBallBoundary<1>:: +// get_normals_at_vertices (const Triangulation<1,1>::face_iterator &, +// Boundary<1,1>::FaceVertexNormals &) const; +// template <> +// Point<1> +// HalfHyperBallBoundary<1>:: +// get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const; +// template <> +// void +// HalfHyperBallBoundary<1>:: +// get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &, +// std::vector > &) const; +// template <> +// void +// HalfHyperBallBoundary<1>:: +// get_normals_at_vertices (const Triangulation<1,1>::face_iterator &, +// Boundary<1,1>::FaceVertexNormals &) const; +// template <> +// Point<1> +// HalfHyperShellBoundary<1>:: +// get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const; +// template <> +// void +// HalfHyperShellBoundary<1>:: +// get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &, +// std::vector > &) const; +// template <> +// void +// HalfHyperShellBoundary<1>:: +// get_normals_at_vertices (const Triangulation<1,1>::face_iterator &, +// Boundary<1,1>::FaceVertexNormals &) const; +// +// +// #endif // DOXYGEN +// +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/include/deal.II/grid/tria.h b/deal.II/include/deal.II/grid/tria.h index adb735e58a..db9f3669d6 100644 --- a/deal.II/include/deal.II/grid/tria.h +++ b/deal.II/include/deal.II/grid/tria.h @@ -42,10 +42,13 @@ DEAL_II_NAMESPACE_OPEN template class Boundary; template class StraightBoundary; +template class Manifold; +template class FlatManifold; template class TriaAccessor; template class TriaAccessor<0,1,spacedim>; + namespace internal { namespace Triangulation @@ -1242,6 +1245,15 @@ public: */ static const StraightBoundary straight_boundary; + /** + * Default manifold object. This is used + * for those objects for which no + * manifold description has been explicitly + * set using set_manifold(). + */ + static const FlatManifold flat_manifold; + + /** * Declare some symbolic names * for mesh smoothing @@ -1808,7 +1820,13 @@ public: void set_boundary (const types::manifold_id number, const Boundary &boundary_object); - /** + + void set_manifold (const types::manifold_id number, + const Manifold &manifold_object); + + + +/** * Reset those parts of the boundary with * the given number to use a straight * boundary approximation. This is the @@ -1823,6 +1841,8 @@ public: */ void set_boundary (const types::manifold_id number); + void set_manifold (const types::manifold_id number); + /** * Return a constant reference to * a boundary object used for @@ -1836,6 +1856,9 @@ public: */ const Boundary &get_boundary (const types::manifold_id number) const; + + const Manifold &get_manifold (const types::manifold_id number) const; + /** * Returns a vector containing * all boundary indicators @@ -3579,6 +3602,16 @@ private: */ std::map , Triangulation > > boundary; + + /** + * Collection of manifold + * objects. We store only + * objects, which are not + * of type StraightBoundary. + */ + std::map , Triangulation > > manifold; + + /** * Flag indicating whether * anisotropic refinement took diff --git a/deal.II/include/deal.II/grid/tria_accessor.h b/deal.II/include/deal.II/grid/tria_accessor.h index 984374f98f..6c2dcc55b9 100644 --- a/deal.II/include/deal.II/grid/tria_accessor.h +++ b/deal.II/include/deal.II/grid/tria_accessor.h @@ -37,6 +37,7 @@ template class TriaIterator; template class TriaActiveIterator; template class Boundary; +template class Manifold; namespace internal @@ -1161,6 +1162,18 @@ public: * for the boundary object. */ const Boundary &get_boundary () const; + + /** + * Return a constant reference to a + * manifold object used for this + * object. This function is a shortcut to + * retrieving the manifold indicator + * using manifold_indicator() and then + * asking the + * Triangulation::get_manifold() function + * for the manifold object. + */ + const Manifold &get_manifold () const; /** * @} 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 6e153eb075..e7d0d57325 100644 --- a/deal.II/include/deal.II/grid/tria_accessor.templates.h +++ b/deal.II/include/deal.II/grid/tria_accessor.templates.h @@ -1959,6 +1959,20 @@ TriaAccessor::get_boundary () const } + +template +const Manifold & +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]; + + return this->tria->get_manifold(mi); +} + + template types::manifold_id TriaAccessor::manifold_id () const diff --git a/deal.II/source/grid/CMakeLists.txt b/deal.II/source/grid/CMakeLists.txt index 5e13fa938c..4968dd0fe6 100644 --- a/deal.II/source/grid/CMakeLists.txt +++ b/deal.II/source/grid/CMakeLists.txt @@ -28,6 +28,8 @@ SET(_src grid_reordering.cc grid_tools.cc intergrid_map.cc + manifold.cc + manifold_lib.cc persistent_tria.cc tria_accessor.cc tria_boundary.cc @@ -45,6 +47,7 @@ SET(_inst grid_refinement.inst.in grid_tools.inst.in intergrid_map.inst.in + manifold.inst.in tria_accessor.inst.in tria_boundary.inst.in tria_boundary_lib.inst.in diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc new file mode 100644 index 0000000000..5135baf965 --- /dev/null +++ b/deal.II/source/grid/manifold.cc @@ -0,0 +1,72 @@ +// --------------------------------------------------------------------- +// $Id: tria_boundary.cc 30130 2013-07-23 13:01:18Z heltai $ +// +// Copyright (C) 1998 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + + +/* -------------------------- Manifold --------------------- */ + + +template +Manifold::~Manifold () +{} + + + + +template +Point +Manifold::project_to_manifold (const Point &candidate) const +{ + Assert (false, ExcPureFunctionCalled()); + return Point(); +} + + +/* -------------------------- FlatManifold --------------------- */ + + +template +FlatManifold::FlatManifold () +{} + + +template +Point +FlatManifold:: +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())); + + Point p; + for(unsigned int i=0; i; + template class FlatManifold; + } + + + diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc index a29e247a5e..9fc1ae83ed 100644 --- a/deal.II/source/grid/tria.cc +++ b/deal.II/source/grid/tria.cc @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -43,6 +44,93 @@ DEAL_II_NAMESPACE_OPEN +/** We create an anonymous namespace to make a helper function which + would return default points and weights contained in lines, quads + and hexes which are used to query the manifold objects. + */ +namespace { + + void get_default_points_and_weights(std::vector > &sp, + std::vector &wp, + const CellAccessor<3,3>& obj) + { + const int dim = 3; + + int np = GeometryInfo::vertices_per_cell+ + GeometryInfo::lines_per_cell+ + GeometryInfo::faces_per_cell; + sp.resize(np); + wp.resize(np); + unsigned int j=0; + for(unsigned int i=0; i::vertices_per_cell; ++i, ++j) + { + sp[j] = obj.vertex(i); + wp[j] = 1.0/128.0; + } + for(unsigned int i=0; i::lines_per_cell; ++i, ++j) + { + sp[j] = obj.line(i)->child(0)->vertex(1); + wp[j] = 7.0/192.0; + } + for(unsigned int i=0; i::faces_per_cell; ++i, ++j) + { + sp[j] = obj.face(i)->isotropic_child(0)->vertex(3); + wp[j] = 1.0/12.0; + } + } + + + template + void get_default_points_and_weights(std::vector > &sp, + std::vector &wp, + const OBJECT& obj, bool with_laplace = false) + { + const int dim = OBJECT::AccessorType::structure_dimension; + AssertDimension(spacedim, OBJECT::AccessorType::space_dimension); + switch(dim) + { + case 1: + sp.resize(2); + wp.resize(2); + sp[0] = obj->vertex(0); wp[0] = .5; + sp[1] = obj->vertex(1); wp[0] = .5; + break; + case 2: + sp.resize(8); + wp.resize(8); + sp[0] = obj->vertex(0); + sp[1] = obj->vertex(1); + sp[2] = obj->vertex(2); + sp[3] = obj->vertex(3); + + sp[4] = obj->line(0)->has_children() ? + obj->line(0)->child(0)->vertex(1) : + obj->line(0)->center(); + sp[5] = obj->line(1)->has_children() ? + obj->line(1)->child(0)->vertex(1) : + obj->line(1)->center(); + sp[6] = obj->line(2)->has_children() ? + obj->line(2)->child(0)->vertex(1) : + obj->line(2)->center(); + sp[7] = obj->line(3)->has_children() ? + obj->line(3)->child(0)->vertex(1) : + obj->line(3)->center(); + if(with_laplace) + { + std::fill(wp.begin(), wp.begin()+4, 1.0/16.0); + std::fill(wp.begin()+4, wp.end(), 3.0/16.0); + } + else + std::fill(wp.begin(), wp.end(), 1.0/8.0); + break; + default: + Assert(false, ExcInternalError()); + break; + } + } +} + + template CellData::CellData() { @@ -3960,8 +4048,7 @@ namespace internal // not the case, then // we need to ask a // boundary object - if (dim == spacedim) - { + { // in a first // step, compute // the location @@ -3970,12 +4057,12 @@ namespace internal // average of the // 8 surrounding // vertices - Point new_point; - for (unsigned int i=0; i<8; ++i) - new_point += triangulation.vertices[new_vertices[i]]; - new_point /= 8.0; - - triangulation.vertices[next_unused_vertex] = new_point; + std::vector > ps(8); + std::vector ws(8, 1.0/8); + for (unsigned int i=0; i<8; ++i) + ps[i] = triangulation.vertices[new_vertices[i]]; + triangulation.vertices[next_unused_vertex] = + cell->get_manifold().get_new_point(ps, ws); // if the user_flag is set, i.e. if the // cell is at the boundary, use a @@ -4014,52 +4101,26 @@ namespace internal boundary_face=GeometryInfo::faces_per_cell+1; } - if (boundary_face::faces_per_cell) - // reset the cell's middle vertex - // to the middle of the straight - // connection between the new - // points on this face and on the - // opposite face - triangulation.vertices[next_unused_vertex] - = 0.5*(cell->face(boundary_face) - ->child(0)->vertex(1)+ - cell->face(GeometryInfo + if (boundary_face::faces_per_cell) + { + // reset the cell's middle vertex + // to the middle of the straight + // connection between the new + // points on this face and on the + // opposite face + std::vector > sp(2); + std::vector ws(2, .5); + sp[0] = cell->face(boundary_face) + ->child(0)->vertex(1); + sp[1] = cell->face(GeometryInfo ::opposite_face[boundary_face]) - ->child(0)->vertex(1)); - } - } - else - { - // if this quad - // lives in a - // higher - // dimensional - // space then we - // don't need to - // worry if it is - // at the - // boundary of - // the manifold - // -- we always - // have to use - // the boundary - // object anyway; - // so ignore - // whether the - // user flag is - // set or not - cell->clear_user_flag(); - - - // new vertex is - // placed on the - // surface according - // to the information - // stored in the - // boundary class - triangulation.vertices[next_unused_vertex] = - cell->get_boundary().get_new_point_on_quad (cell); - } + ->child(0)->vertex(1); + + triangulation.vertices[next_unused_vertex] + = cell->get_manifold().get_new_point(ps, ws); + } + } + } } @@ -4316,6 +4377,9 @@ namespace internal const bool /*check_for_distorted_cells*/) { const unsigned int dim = 1; + + std::vector > sp(2); + std::vector wp(2); // check whether a new level is // needed we have to check for this @@ -4444,12 +4508,10 @@ namespace internal // midpoint; otherwise we // have to ask the manifold // object - if (dim == spacedim) - triangulation.vertices[next_unused_vertex] = - (cell->vertex(0) + cell->vertex(1)) / 2; - else - triangulation.vertices[next_unused_vertex] = - cell->get_boundary().get_new_point_on_line(cell); + get_default_points_and_weights(sp, wp, cell); + + triangulation.vertices[next_unused_vertex] = + cell->get_manifold().get_new_point(sp, wp); triangulation.vertices_used[next_unused_vertex] = true; @@ -4605,6 +4667,15 @@ namespace internal { const unsigned int dim = 2; + // Helper vectors of points + // and weights, to compute + // the locations of the newly + // created vertices (through + // the manifold objects) + std::vector > sp; + std::vector wp; + + // check whether a new level is // needed we have to check for this // on the highest level only (on @@ -4862,8 +4933,9 @@ namespace internal // Ask the manifold // where to put the // new line - triangulation.vertices[next_unused_vertex] - = line->get_boundary().get_new_point_on_line(line); + get_default_points_and_weights(sp, wp, line); + triangulation.vertices[next_unused_vertex] = + line->get_manifold().get_new_point(sp, wp); // if (spacedim == dim) @@ -5043,6 +5115,9 @@ namespace internal { const unsigned int dim = 3; + std::vector > sp; + std::vector wp; + // this function probably // also works for spacedim>3 // but it isn't tested. it @@ -5423,8 +5498,10 @@ namespace internal ExcTooFewVerticesAllocated()); triangulation.vertices_used[next_unused_vertex] = true; - triangulation.vertices[next_unused_vertex] - = line->get_boundary().get_new_point_on_line(line); + get_default_points_and_weights(sp, wp, line); + triangulation.vertices[next_unused_vertex] = + line->get_manifold().get_new_point(sp, wp); + // if (line->at_boundary()) // triangulation.vertices[next_unused_vertex] @@ -6033,9 +6110,10 @@ namespace internal // the domain, so we need not // care about boundary quads // here - triangulation.vertices[next_unused_vertex] - = middle_line->get_boundary().get_new_point_on_line(middle_line); -// = (middle_line->vertex(0) + middle_line->vertex(1)) / 2; + get_default_points_and_weights(sp, wp, middle_line); + triangulation.vertices[next_unused_vertex] = + middle_line->get_manifold().get_new_point(sp, wp); + triangulation.vertices_used[next_unused_vertex] = true; // now search a slot for the two @@ -6099,65 +6177,53 @@ namespace internal // isotropic refinement Assert(quad_ref_case==RefinementCase::no_refinement, ExcInternalError()); - // set the middle vertex - // appropriately - if (quad->at_boundary() || quad->manifold_id()!= numbers::invalid_manifold_id) - triangulation.vertices[next_unused_vertex] - = quad->get_boundary().get_new_point_on_quad(quad); -// = triangulation.get_boundary(quad->boundary_indicator()).get_new_point_on_quad (quad); - else - // it might be that the - // quad itself is not - // at the boundary, but - // that one of its lines - // actually is. in this - // case, the newly - // created vertices at - // the centers of the - // lines are not - // necessarily the mean - // values of the - // adjacent vertices, - // so do not compute - // the new vertex as - // the mean value of - // the 4 vertices of - // the face, but rather - // as a weighted mean - // value of the 8 - // vertices which we - // already have (the - // four old ones, and - // the four ones - // inserted as middle - // points for the four - // lines). summing up - // some more points is - // generally cheaper - // than first asking - // whether one of the - // lines is at the - // boundary - // - // note that the exact - // weights are chosen - // such as to minimize - // the distortion of - // the four new quads - // from the optimal - // shape; their - // derivation and - // values is copied - // over from the - // @p{MappingQ::set_laplace_on_vector} - // function - triangulation.vertices[next_unused_vertex] - = (quad->vertex(0) + quad->vertex(1) + - quad->vertex(2) + quad->vertex(3) + - 3*(quad->line(0)->child(0)->vertex(1) + - quad->line(1)->child(0)->vertex(1) + - quad->line(2)->child(0)->vertex(1) + - quad->line(3)->child(0)->vertex(1)) ) / 16; + // // set the middle vertex + // // appropriately + // get_default_points_and_weights(sp, wp, middle_line); + // triangulation.vertices[next_unused_vertex] = + // middle_line->get_manifold().get_new_point(sp, wp); + + if (quad->at_boundary() || quad->manifold_id()!= numbers::invalid_manifold_id) + { + get_default_points_and_weights(sp, wp, quad); + triangulation.vertices[next_unused_vertex] + = quad->get_manifold().get_new_point(sp, wp); + } + else + { + + // it might be that the quad itself is not at + // the boundary, but that one of its lines + // actually is. in this case, the newly created + // vertices at the centers of the lines are not + // necessarily the mean values of the adjacent + // vertices, so do not compute the new vertex as + // the mean value of the 4 vertices of the face, + // but rather as a weighted mean value of the 8 + // vertices which we already have (the four old + // ones, and the four ones inserted as middle + // points for the four lines). summing up some + // more points is generally cheaper than first + // asking whether one of the lines is at the + // boundary + // + // note that the exact weights are chosen such + // as to minimize the distortion of the four new + // quads from the optimal shape; their + // derivation and values is copied over from the + // @p{MappingQ::set_laplace_on_vector} function + get_default_points_and_weights(sp, wp, quad, true); + triangulation.vertices[next_unused_vertex] + = quad->get_manifold().get_new_point(sp, wp); + + // triangulation.vertices[next_unused_vertex] + // = (quad->vertex(0) + quad->vertex(1) + + // quad->vertex(2) + quad->vertex(3) + + // 3*(quad->line(0)->child(0)->vertex(1) + + // quad->line(1)->child(0)->vertex(1) + + // quad->line(2)->child(0)->vertex(1) + + // quad->line(3)->child(0)->vertex(1)) ) / 16; + } triangulation.vertices_used[next_unused_vertex] = true; @@ -8423,26 +8489,29 @@ namespace internal // vertex at the center of // the quad is weighted (see // above) - triangulation.vertices[next_unused_vertex] = Point(); - // first add corners of hex - for (unsigned int vertex=0; - vertex::vertices_per_cell; ++vertex) - triangulation.vertices[next_unused_vertex] += hex->vertex(vertex) / 128; - // now add center of lines - for (unsigned int line=0; - line::lines_per_cell; ++line) - triangulation.vertices[next_unused_vertex] += hex->line(line)->child(0)->vertex(1) * - 7./192.; - // finally add centers of - // faces. note that vertex 3 - // of child 0 is an invariant - // with respect to the face - // orientation, flip and - // rotation - for (unsigned int face=0; - face::faces_per_cell; ++face) - triangulation.vertices[next_unused_vertex] += hex->face(face)->isotropic_child(0)->vertex(3) * - 1./12.; + get_default_points_and_weights(sp, wp, hex); + triangulation.vertices[next_unused_vertex] = + hex->get_manifold().get_new_point(sp, wp); + + // // first add corners of hex + // for (unsigned int vertex=0; + // vertex::vertices_per_cell; ++vertex) + // triangulation.vertices[next_unused_vertex] += hex->vertex(vertex) / 128; + // // now add center of lines + // for (unsigned int line=0; + // line::lines_per_cell; ++line) + // triangulation.vertices[next_unused_vertex] += hex->line(line)->child(0)->vertex(1) * + // 7./192.; + // // finally add centers of + // // faces. note that vertex 3 + // // of child 0 is an invariant + // // with respect to the face + // // orientation, flip and + // // rotation + // for (unsigned int face=0; + // face::faces_per_cell; ++face) + // triangulation.vertices[next_unused_vertex] += hex->face(face)->isotropic_child(0)->vertex(3) * + // 1./12.; // set the data of the // six lines. first @@ -9501,6 +9570,10 @@ const StraightBoundary Triangulation::straight_boundary = StraightBoundary(); +template +const FlatManifold +Triangulation::flat_manifold = FlatManifold(); + template const unsigned int @@ -9610,6 +9683,18 @@ Triangulation::set_boundary (const types::manifold_id m_number, +template +void +Triangulation::set_manifold (const types::manifold_id m_number, + const Manifold &manifold_object) +{ + Assert(m_number < numbers::invalid_manifold_id, + ExcIndexRange(m_number,0,numbers::invalid_manifold_id)); + + manifold[m_number] = &manifold_object; +} + + template void Triangulation::set_boundary (const types::manifold_id m_number) @@ -9622,6 +9707,18 @@ Triangulation::set_boundary (const types::manifold_id m_number) } +template +void +Triangulation::set_manifold (const types::manifold_id m_number) +{ + Assert(m_number < numbers::invalid_manifold_id, + ExcIndexRange(m_number,0,numbers::invalid_manifold_id)); + + //delete the entry located at number. + manifold.erase(m_number); +} + + template const Boundary & @@ -9650,6 +9747,35 @@ Triangulation::get_boundary (const types::manifold_id m_number) c } +template +const Manifold & +Triangulation::get_manifold (const types::manifold_id m_number) const +{ + Assert(m_number < numbers::invalid_manifold_id, + ExcIndexRange(m_number,0,numbers::invalid_manifold_id)); + + //look, if there is a boundary stored at + //boundary_id number. + typename std::map, Triangulation > >::const_iterator it + = manifold.find(m_number); + + if (it != manifold.end()) + { + //if we have found an entry, return it + return *(it->second); + } + else + { + //if we have not found an entry + //connected with number, we return + //flat_manifold + return flat_manifold; + } +} + + + + template std::vector Triangulation::get_boundary_indicators () const diff --git a/deal.II/source/grid/tria.inst.in b/deal.II/source/grid/tria.inst.in index 66e61c2864..c1f3faa103 100644 --- a/deal.II/source/grid/tria.inst.in +++ b/deal.II/source/grid/tria.inst.in @@ -22,6 +22,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS template class Triangulation; #endif #if deal_II_dimension == deal_II_space_dimension - template class CellData; + template struct CellData; #endif } -- 2.39.5