--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <deal.II/base/config.h>
+#include <deal.II/base/subscriptor.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/base/subscriptor.h>
+#include <deal.II/base/point.h>
+
+#include <deal.II/grid/tria.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <int dim, int space_dim> 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<dim> new_vertex;
+ *
+ * for(unsigned int i=0; i<surrounding_points.size(); ++i)
+ * new_vertex += weights[i]*surrounding_points[i];
+ * ...
+ * @endcode
+ *
+ * Each different manifold description will override the function above,
+ * by assuming that the surrounding points are on the manifold, and computing
+ * whatever is meaninful out of the weighted average of the surrounding_points.
+ *
+ * @ingroup manifold
+ * @author Luca Heltai, 2013
+ */
+template <int spacedim>
+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<spacedim>
+ get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &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<spacedim> project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
+ const Point<spacedim> &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
+ * <tt>quad->line(i)->child(j)</tt>, <tt>i=0...3</tt>, <tt>j=0,1</tt>.
+ *
+ * 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<spacedim>, which represent the
+ * Euclidean space of dimension #spacedim.
+ *
+ * @ingroup manifold
+ *
+ * @author Luca Heltai, 2013
+ */
+template <int spacedim>
+class FlatManifold: public Manifold<spacedim>
+{
+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<spacedim> periodicity=Point<spacedim>());
+
+ /**
+ * 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<spacedim>
+ get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &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<spacedim> project_to_manifold (const std::vector<Point<spacedim> > &points,
+ const Point<spacedim> &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<spacedim> periodicity;
+};
+
+
+/**
+ * A chart of dimension dim, which is part of a Manifold<spacedim>.
+ * 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 <int dim, int spacedim>
+class ManifoldChart: public Manifold<spacedim>
+{
+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<dim> periodicity=Point<dim>());
+
+ /**
+ * 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<spacedim>
+ get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &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<dim>
+ pull_back(const Point<spacedim> &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<spacedim>
+ push_forward(const Point<dim> &chart_point) const = 0;
+
+ private:
+ /**
+ * The sub_manifold object is used to compute the average of the
+ * points in the chart coordinates system.
+ */
+ FlatManifold<dim> sub_manifold;
+};
+
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
+#ifndef DOXYGEN
+
+#endif // DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <deal.II/base/config.h>
+#include <deal.II/grid/manifold.h>
+
+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 <int spacedim>
+class PolarManifold : public FlatManifold<spacedim>
+{
+ public:
+ /**
+ * The Constructor takes the
+ * center of the polar
+ * coordinates system.
+ */
+ PolarManifold(const Point<spacedim> center = Point<spacedim>());
+
+ /**
+ * 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<spacedim>
+ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const;
+
+
+ /**
+ * The center of the polar
+ * coordinate system.
+ */
+ const Point<spacedim> center;
+};
+
+/**
+ * Manifold description for a cylindrical space coordinate system. In
+ * three dimensions, points are projected on a cylindrical tube along
+ * the <tt>x-</tt>, <tt>y-</tt> or <tt>z</tt>-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 <int spacedim>
+class CylinderManifold : public FlatManifold<spacedim>
+{
+ public:
+ /**
+ * Constructor. Using default
+ * values for the constructor
+ * arguments yields a cylindrical
+ * manifold tube along the x-axis
+ * (<tt>axis=0</tt>). Choose
+ * <tt>axis=1</tt> or
+ * <tt>axis=2</tt> 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<spacedim> &direction,
+ const Point<spacedim> &point_on_axis);
+
+ /**
+ * Refer to the general
+ * documentation of this class
+ * and the documentation of the
+ * base class.
+ */
+ virtual Point<spacedim>
+ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const;
+
+
+ private:
+ /**
+ * Direction of the cylindrical
+ * coordinate system.
+ */
+ const Point<spacedim> direction;
+
+ /**
+ * One point on the axis of the
+ * cylindrical coordinate systems
+ */
+ const Point<spacedim> point_on_axis;
+
+ /**
+ * Returns the axis vector
+ * associated with the given axis
+ * number.
+ */
+ static Point<spacedim> get_axis_vector (const unsigned int axis);
+};
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
#include <deal.II/base/thread_management.h>
#include <deal.II/base/point.h>
#include <deal.II/grid/tria.h>
+#include <deal.II/grid/manifold.h>
DEAL_II_NAMESPACE_OPEN
template <int dim, int space_dim> class Triangulation;
+template <int spacedim> class Manifold;
+
/**
* @author Wolfgang Bangerth, 1999, 2001, 2009, Ralf Hartmann, 2001, 2008
*/
template <int dim, int spacedim=dim>
-class Boundary : public Subscriptor
+class Boundary : public Manifold<spacedim>
{
public:
grid_reordering.cc
grid_tools.cc
intergrid_map.cc
+ manifold.cc
+ manifold_lib.cc
persistent_tria.cc
tria_accessor.cc
tria_boundary.cc
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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <deal.II/base/point.h>
+#include <deal.II/grid/manifold.h>
+#include <cmath>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+/* -------------------------- Manifold --------------------- */
+
+
+template <int spacedim>
+Manifold<spacedim>::~Manifold ()
+{}
+
+
+template <int spacedim>
+Point<spacedim>
+Manifold<spacedim>::project_to_manifold (const std::vector<Point<spacedim> > &,
+ const Point<spacedim> &candidate) const
+{
+ Assert (false, ExcPureFunctionCalled());
+ return candidate;
+}
+
+/* -------------------------- FlatManifold --------------------- */
+
+
+template <int spacedim>
+FlatManifold<spacedim>::FlatManifold (const Point<spacedim> periodicity) :
+ periodicity(periodicity)
+{}
+
+template <int spacedim>
+Point<spacedim>
+FlatManifold<spacedim>::
+get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const
+{
+ AssertDimension(surrounding_points.size(), weights.size());
+
+#ifdef DEBUG
+ double sum=0;
+ for(unsigned int i=0; i<weights.size(); ++i)
+ sum+= weights[i];
+ Assert(std::abs(sum-1.0) < 1e-10, ExcMessage("Weights should sum to 1!"));
+#endif
+
+
+ Point<spacedim> p;
+ Point<spacedim> dp;
+ bool check_period = (periodicity.norm() != 0);
+ for(unsigned int i=0; i<surrounding_points.size(); ++i) {
+ if(check_period) {
+ for(unsigned int d=0; d<spacedim; ++d)
+ if(periodicity[d] > 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<spacedim; ++d)
+ if(periodicity[d] > 0)
+ p[d] = std::fmod(p[d], periodicity[d]);
+
+ return project_to_manifold(surrounding_points, p);
+}
+
+template <int spacedim>
+Point<spacedim>
+FlatManifold<spacedim>::project_to_manifold (const std::vector<Point<spacedim> > & vertices,
+ const Point<spacedim> &candidate) const
+{
+ return candidate;
+}
+
+
+/* -------------------------- ManifoldChart --------------------- */
+
+template <int dim, int spacedim>
+ManifoldChart<dim,spacedim>::~ManifoldChart ()
+{}
+
+template <int dim, int spacedim>
+ManifoldChart<dim,spacedim>::ManifoldChart (const Point<dim> periodicity):
+ sub_manifold(periodicity)
+{}
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+ManifoldChart<dim, spacedim>::
+get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const
+{
+ AssertDimension(surrounding_points.size(), weights.size());
+ std::vector<Point<dim> > chart_points(surrounding_points.size());
+
+ for(unsigned int i=0; i<surrounding_points.size(); ++i)
+ chart_points[i] = pull_back(surrounding_points[i]);
+
+ Point<dim> p_chart = sub_manifold.get_new_point(chart_points, weights);
+
+ return push_forward(p_chart);
+}
+
+
+// explicit instantiations
+#include "manifold.inst"
+
+DEAL_II_NAMESPACE_CLOSE
+
--- /dev/null
+// ---------------------------------------------------------------------
+// $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<deal_II_space_dimension>;
+ template class FlatManifold<deal_II_space_dimension>;
+ template class ManifoldChart<deal_II_space_dimension,deal_II_space_dimension>;
+#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
+
+ }
+
+
+
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/base/tensor.h>
+#include <cmath>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <int spacedim>
+PolarManifold<spacedim>::PolarManifold(const Point<spacedim> center):
+ center(center)
+{}
+
+
+
+template <int spacedim>
+Point<spacedim>
+PolarManifold<spacedim>::
+get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const
+{
+ Assert(surrounding_points.size() == weights.size(),
+ ExcDimensionMismatch(surrounding_points.size(), weights.size()));
+
+ double radius = 0;
+ for(unsigned int i=0; i<surrounding_points.size(); ++i)
+ radius += weights[i]*(surrounding_points[i]-center).norm();
+
+ Point<spacedim> p = FlatManifold<spacedim>::get_new_point(surrounding_points, weights);
+ p = p-center;
+ p = p/p.norm()*radius+center;
+ return p;
+}
+
+// ============================================================
+// CylinderManifold
+// ============================================================
+
+template <int spacedim>
+CylinderManifold<spacedim>::CylinderManifold(const int axis):
+ direction (get_axis_vector (axis)),
+ point_on_axis (Point<spacedim>())
+{}
+
+template <int spacedim>
+CylinderManifold<spacedim>::CylinderManifold(const Point<spacedim> &direction,
+ const Point<spacedim> &point_on_axis):
+ direction (direction),
+ point_on_axis (point_on_axis)
+{}
+
+
+
+template <int spacedim>
+Point<spacedim>
+CylinderManifold<spacedim>::get_axis_vector (const unsigned int axis)
+{
+
+ Assert (axis < spacedim, ExcIndexRange (axis, 0, spacedim));
+
+ Point<spacedim> axis_vector;
+
+ axis_vector[axis] = 1;
+
+ return axis_vector;
+}
+
+
+template <int spacedim>
+Point<spacedim>
+CylinderManifold<spacedim>::
+get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const
+{
+ Assert(surrounding_points.size() == weights.size(),
+ ExcDimensionMismatch(surrounding_points.size(), weights.size()));
+
+ // compute a proposed new point
+ Point<spacedim> middle = FlatManifold<spacedim>::get_new_point(surrounding_points, weights);
+
+ double radius = 0;
+ Point<spacedim> on_plane;
+
+ for(unsigned int i=0; i<surrounding_points.size(); ++i)
+ {
+ on_plane = surrounding_points[i]-point_on_axis;
+ on_plane = on_plane - (on_plane*direction) * direction;
+ radius += weights[i]*on_plane.norm();
+ }
+
+ // we then have to project this
+ // point out to the given radius
+ // from the axis. to this end, we
+ // have to take into account the
+ // offset point_on_axis and the
+ // direction of the axis
+ const Point<spacedim> 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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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<deal_II_dimension>;
+ template class CylinderManifold<deal_II_dimension>;
+// template class ConeBoundary<deal_II_dimension>;
+ }
+
+
projects / dealii-svn.git / commitdiff