From: heltai Date: Thu, 17 Jul 2014 17:53:51 +0000 (+0000) Subject: Added Manifold class, from which Boundary is derived. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9a55d4b85c5b5cfe5b40f8745b47e8accebabce5;p=dealii-svn.git Added Manifold class, from which Boundary is derived. git-svn-id: https://svn.dealii.org/branches/branch_boundary_to_manifold@33193 0785d39b-7218-0410-832d-ea1e28bc413d --- 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..063a7c7cb4 --- /dev/null +++ b/deal.II/include/deal.II/grid/manifold.h @@ -0,0 +1,356 @@ +// --------------------------------------------------------------------- +// $Id: tria_boundary.h 30130 2013-07-23 13:01:18Z heltai $ +// +// Copyright (C) 2013, 2014 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 + + +/*---------------------------- manifold.h ---------------------------*/ + +#include +#include +#include +#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; + + /** + * Given a point which lies close to the given manifold, it modifies + * it and projects it to manifold itself. 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 manfold remain on the manifold. + * + * Derived classes do not need to implement this function unless + * mesh smoothing algorithms are used with a particular boundary + * object. + */ + virtual + Point project_to_manifold (const std::vector > &surrounding_points, + const Point &candidate) const; + + + virtual + Point<2> + get_new_point_on_line (const typename Triangulation<1,2>::line_iterator &line) const; + + /** + * Return the point which shall become the new middle vertex of the two + * children of a regular line. In 2D, this line is a line at the boundary, + * while in 3d, it is bounding a face at the boundary (the lines therefore + * is also on the boundary). + */ + virtual + Point<2> + get_new_point_on_line (const typename Triangulation<2,2>::line_iterator &line) const; + + virtual + Point<3> + get_new_point_on_line (const typename Triangulation<2,3>::line_iterator &line) const; + + virtual + Point<3> + get_new_point_on_line (const typename Triangulation<3,3>::line_iterator &line) const; + + /** + * Return the point which shall become the common point of the four children + * of a quad at the boundary in three or more spatial dimensions. This + * function therefore is only useful in at least three dimensions and should + * not be called for lower dimensions. + * + * This function is called after the four lines bounding the given @p quad + * are refined, so you may want to use the information provided by + * quad->line(i)->child(j), i=0...3, j=0,1. + * + * Because in 2D, this function is not needed, it is not made pure virtual, + * to avoid the need to overload it. The default implementation throws an + * error in any case, however. + */ + virtual + Point<2> + get_new_point_on_quad (const typename Triangulation<2,2>::quad_iterator &quad) const; + + virtual + Point<3> + get_new_point_on_quad (const typename Triangulation<2,3>::quad_iterator &quad) const; + + virtual + Point<3> + get_new_point_on_quad (const typename Triangulation<3,3>::quad_iterator &quad) const; + + virtual + Point<3> + get_new_point_on_hex (const typename Triangulation<3,3>::hex_iterator &hex) const; + +private: + +}; + + + +/** + * Specialization of Manifold, which represent the + * Euclidean space of dimension #spacedim. + * + * @ingroup manifold + * + * @author Luca Heltai, 2013 + */ +template +class FlatManifold: public Manifold +{ +public: + /** + * Default constructor. Some compilers require this for some + * reasons. The The optional argument can be used to specify the + * periodicity of the spacedim-dimensional manifold (one period per + * direction). A peridicity value of zero means that along that + * direction there is no peridicity. By default no periodicity is + * assumed. + * + * Periodicity affects the way a middle point is computed. It is + * assumed that if two points are more than half period distant, + * then the distance should be computed by crossing the periodicity + * boundary, i.e., the average is computed by adding a full period + * to the sum of the two. For example, if along direction 0 we have + * 2*pi periodicity, then the average of (2*pi-eps) and (eps) is not + * pi, but 2*pi (or zero), since, on a periodic manifold, these two + * points are at distance 2*eps and not (2*pi-eps) + */ + FlatManifold (const Point periodicity=Point()); + + /** + * Let the new point be the average sum of surrounding vertices. + * + * This particular implementation constructs the weighted average of + * the surrounding points, and then calss internally the function + * project_to_manifold. The reason why we do it this way, is to + * allow lazy programmers to implement only the project_to_manifold + * function for their own Manifold classes which are small (or + * trivial) perturbations of a flat manifold. For most simple + * geometries, it is possible to get reasonable results by deriving + * your own Manifold class from FlatManifold, and write a new + * interface only for the project_to_manifold function. + */ + virtual Point + get_new_point(const std::vector > &surrounding_points, + const std::vector &weights) const; + + /** + * Project to FlatManifold: do nothing. Note however that this + * function can be overloaded by derived classes, which will then + * benefit from the logic behind the get_new_point class which are + * often very similar (if not identical) to the one implemented in + * this class. + */ + virtual + Point project_to_manifold (const std::vector > &points, + const Point &candidate) const; + +private: + /** + * The periodicity of this Manifold. Periodicity affects the way a + * middle point is computed. It is assumed that if two points are + * more than half period distant, then the distance should be + * computed by crossing the periodicity boundary, i.e., the average + * is computed by adding a full period to the sum of the two. For + * example, if along direction 0 we have 2*pi periodicity, then the + * average of (2*pi-eps) and (eps) is not pi, but 2*pi (or zero), + * since, on a periodic manifold, these two points are at distance + * 2*eps and not (2*pi-eps). + * + * A periodicity 0 along one direction means no periodicity. This is + * the default value for all directions. + */ + const Point periodicity; +}; + + +/** + * A chart of dimension dim, which is part of a Manifold. + * This object specializes a Manifold of dimension dim embedded in a + * manifold of dimension spacedim, for which you have explicit + * pull_back and push_forward transformations. + * + * This is an helper class which is useful when you have an explicit + * map from an Euclidean space of dimension dim to an Euclidean + * space of dimension spacedim which represent your manifold, i.e., + * when your manifold \f$\mathcal{M}\f$ can be represented by a map + * \f[ + * F: \mathcal{B} \subset R^{\text{dim}} \mapsto \mathcal{M} + * \subset R^{\text{spacedim}} + * \f] + * (the push_forward() function) + * which admits the inverse transformation + * \f[ + * F^{-1}: \mathcal{M} + * \subset R^{\text{spacedim}} \mapsto + * \mathcal{B} \subset R^{\text{dim}} + * \f] + * (the pull_back() function). + * + * The get_new_point() function of the ManifoldChart class is + * implemented by calling the pull_back() method for all + * #surrounding_points, computing their weighted average in the dim + * Euclidean space, and calling the push_forward() method with the + * resulting point, i.e., + * \f[ + * p^{\text{new}} = F(\sum_i w_i F^{-1}(p_i)). + * \f] + * + * Derived classes are required to implement the push_forward() and + * the pull_back() methods. + * + * @ingroup manifold + * + * @author Luca Heltai, 2013 + */ +template +class ManifoldChart: public Manifold +{ +public: + /** + * Constructor. The optional argument can be used to specify the + * periodicity of the dim-dimensional manifold (one period per + * direction). A peridicity value of zero means that along that + * direction there is no peridicity. By default no periodicity is + * assumed. + * + * Periodicity affects the way a middle point is computed. It is + * assumed that if two points are more than half period distant, + * then the distance should be computed by crossing the periodicity + * boundary, i.e., then the average is computed by adding a full + * period to the sum of the two. For example, if along direction 0 + * we have 2*pi periodicity, then the average of (2*pi-eps) and + * (eps) is not pi, but 2*pi (or zero), since, on the manifold, + * these two points are at distance 2*eps and not (2*pi-eps) + */ + ManifoldChart(const Point periodicity=Point()); + + /** + * Destructor. Does nothing here, but needs to be declared to make + * it virtual. + */ + virtual ~ManifoldChart (); + + + /** + * 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; + + /** + * Pull back the given point in spacedim to the Euclidean dim + * dimensional space. + * + * Refer to the general documentation of this class for more + * information. + */ + virtual Point + pull_back(const Point &space_point) const = 0; + + /** + * Given a point in the dim dimensianal Euclidean space, this + * method returns a point on the manifold embedded in the spacedim + * Euclidean space. + * + * Refer to the general documentation of this class for more + * information. + */ + virtual Point + push_forward(const Point &chart_point) const = 0; + + private: + /** + * The sub_manifold object is used to compute the average of the + * points in the chart coordinates system. + */ + FlatManifold sub_manifold; +}; + + + +/* -------------- 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..fe53dc2f15 --- /dev/null +++ b/deal.II/include/deal.II/grid/manifold_lib.h @@ -0,0 +1,165 @@ +// --------------------------------------------------------------------- +// $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 + +/** + * 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; + + + /** + * 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); +}; + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/include/deal.II/grid/tria_boundary.h b/deal.II/include/deal.II/grid/tria_boundary.h index ca948550fc..582e7a63ca 100644 --- a/deal.II/include/deal.II/grid/tria_boundary.h +++ b/deal.II/include/deal.II/grid/tria_boundary.h @@ -26,11 +26,14 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN template class Triangulation; +template class Manifold; + /** @@ -82,7 +85,7 @@ template class Triangulation; * @author Wolfgang Bangerth, 1999, 2001, 2009, Ralf Hartmann, 2001, 2008 */ template -class Boundary : public Subscriptor +class Boundary : public Manifold { public: diff --git a/deal.II/source/grid/CMakeLists.txt b/deal.II/source/grid/CMakeLists.txt index fcac409a1b..fff6739c6d 100644 --- a/deal.II/source/grid/CMakeLists.txt +++ b/deal.II/source/grid/CMakeLists.txt @@ -24,6 +24,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 @@ -41,6 +43,8 @@ SET(_inst grid_refinement.inst.in grid_tools.inst.in intergrid_map.inst.in + manifold.inst.in + manifold_lib.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..4c2c454ada --- /dev/null +++ b/deal.II/source/grid/manifold.cc @@ -0,0 +1,130 @@ +// --------------------------------------------------------------------- +// $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 std::vector > &, + const Point &candidate) const +{ + Assert (false, ExcPureFunctionCalled()); + return candidate; +} + +/* -------------------------- FlatManifold --------------------- */ + + +template +FlatManifold::FlatManifold (const Point periodicity) : + periodicity(periodicity) +{} + +template +Point +FlatManifold:: +get_new_point (const std::vector > &surrounding_points, + const std::vector &weights) const +{ + AssertDimension(surrounding_points.size(), weights.size()); + +#ifdef DEBUG + double sum=0; + for(unsigned int i=0; i p; + Point dp; + bool check_period = (periodicity.norm() != 0); + for(unsigned int i=0; i 0) + dp[d] = std::abs(surrounding_points[i][d]- + surrounding_points[0][d]) > periodicity[d]/2.0 ? + periodicity[d] : 0.0; + } + p += (surrounding_points[i]+dp)*weights[i]; + } + if(check_period) + for(unsigned int d=0; d 0) + p[d] = std::fmod(p[d], periodicity[d]); + + return project_to_manifold(surrounding_points, p); +} + +template +Point +FlatManifold::project_to_manifold (const std::vector > & vertices, + const Point &candidate) const +{ + return candidate; +} + + +/* -------------------------- ManifoldChart --------------------- */ + +template +ManifoldChart::~ManifoldChart () +{} + +template +ManifoldChart::ManifoldChart (const Point periodicity): + sub_manifold(periodicity) +{} + + +template +Point +ManifoldChart:: +get_new_point (const std::vector > &surrounding_points, + const std::vector &weights) const +{ + AssertDimension(surrounding_points.size(), weights.size()); + std::vector > chart_points(surrounding_points.size()); + + for(unsigned int i=0; i p_chart = sub_manifold.get_new_point(chart_points, weights); + + return push_forward(p_chart); +} + + +// explicit instantiations +#include "manifold.inst" + +DEAL_II_NAMESPACE_CLOSE + diff --git a/deal.II/source/grid/manifold.inst.in b/deal.II/source/grid/manifold.inst.in new file mode 100644 index 0000000000..7f99d819b1 --- /dev/null +++ b/deal.II/source/grid/manifold.inst.in @@ -0,0 +1,34 @@ +// --------------------------------------------------------------------- +// $Id: tria_boundary.inst.in 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. +// +// --------------------------------------------------------------------- + + + +for (deal_II_space_dimension : SPACE_DIMENSIONS) + { + template class Manifold; + template class FlatManifold; + template class ManifoldChart; +#if deal_II_space_dimension > 1 + template class ManifoldChart<1,deal_II_space_dimension>; +#endif +#if deal_II_space_dimension > 2 + template class ManifoldChart<2,deal_II_space_dimension>; +#endif + + } + + + diff --git a/deal.II/source/grid/manifold_lib.cc b/deal.II/source/grid/manifold_lib.cc new file mode 100644 index 0000000000..acac37057c --- /dev/null +++ b/deal.II/source/grid/manifold_lib.cc @@ -0,0 +1,134 @@ +// --------------------------------------------------------------------- +// $Id: manifold_lib.cc 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. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +template +PolarManifold::PolarManifold(const Point center): + center(center) +{} + + + +template +Point +PolarManifold:: +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())); + + double radius = 0; + for(unsigned int i=0; i p = FlatManifold::get_new_point(surrounding_points, weights); + p = p-center; + p = p/p.norm()*radius+center; + 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" + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/grid/manifold_lib.inst.in b/deal.II/source/grid/manifold_lib.inst.in new file mode 100644 index 0000000000..11bf4195a9 --- /dev/null +++ b/deal.II/source/grid/manifold_lib.inst.in @@ -0,0 +1,25 @@ +// --------------------------------------------------------------------- +// $Id: tria_boundary_lib.inst.in 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. +// +// --------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class PolarManifold; + template class CylinderManifold; +// template class ConeBoundary; + } + +