Added specific instantiations for invalid accessors. Moved get_default_quadrature into an anonymous namespace.
Made Boundary derived from FlatManifold, and made sure that nothing broke. Added get_new_point_on_cell/hex. Made get_new_point_on_quad not throw an exception in 2d.
Created Manifolds namespace, and made specialization of get_new_point_on_hex only in 3d, instead of multiple specializations for when dim not equal 3.
Added test flat_manifold_01
Fixed 2d problems.
Fixed 3d. Now all tests work, except from machine precision errors.
Added a new test for flat manifold checks on faces.
Completed documentation of manifold and tria_boundary.
Removed commented out functions from Boundary<dim,spacedim>, whose default has been moved to Manifold<dim,spacedim>
Fixed periodicity of FlatManifold, and added a test for it.
Fixed comment on periodicity.
Removed commented instantiations.
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__tria_manifold_h
+#define __deal2__tria_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/point.h>
+#include <deal.II/grid/tria.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <int dim, int space_dim> class Triangulation;
+
+
+/** We collect here some helper functions used in the
+ Manifold<dim,spacedim> classes.
+ */
+namespace Manifolds {
+ /** Given a hex iterator, construct a quadrature with the Laplace
+ weigths, and all relevant points of the hex: vertices, line
+ centers and face centers, which can be called when creating
+ middle vertices in the manifold routines.*/
+ void get_default_quadrature(const TriaIterator<CellAccessor<3, 3> >& hex,
+ Quadrature<3> &quad);
+
+ /** Given a general mesh iterator, construct a quadrature with the
+ Laplace weigths or with uniform weights according the parameter
+ @p with_laplace, and with all relevant points of the iterator:
+ vertices, line centers and/or face centers, which can be called
+ when creating new vertices in the manifold routines.*/
+ template <typename OBJECT, int spacedim>
+ void get_default_quadrature(const OBJECT& obj,
+ Quadrature<spacedim> &quad,
+ bool with_laplace = false);
+}
+
+
+/**
+ * This class is used to represent a manifold to a triangulation.
+ * When a triangulation creates a new vertex on this manifold, it
+ * determines the new vertex' coordinates through the following
+ * function:
+ *
+ * @code
+ * ...
+ * Point<spacedim> new_vertex = manifold.get_new_point (quadrature);
+ * ...
+ * @endcode
+ * @p quadrature is a Quadrature<spacedim> object, which contains a
+ * collection of points in @p spacedim dimension, and a collection of
+ * weights.
+ *
+ * Internally, the get_new_point() function calls the project_to_manifold()
+ * function. This allows end users to only overload project_to_manifold().
+ *
+ * Should a finer control be necessary, then get_new_point() can be
+ * overloaded. For backward compatibility, this function also
+ * offers an interface which is compatible with
+ * Boundary<dim,spacedim>, which are all derived from
+ * FlatManifold<dim,spacedim>, allowing old user codes to keep using
+ * their boundary descriptors as Manifold<dim,spacedim> objects.
+ *
+ * FlatManifold is the specialization from which StraigthBoundary is
+ * derived, where the project_to_manifold() function is the identity.
+ *
+ * @ingroup manifold
+ * @author Luca Heltai, 2014
+ */
+template <int dim, int spacedim=dim>
+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 which make up the quadrature. We use a
+ * quadrature object, which should be filled with the surrounding
+ * points together with appropriate weights.
+ *
+ * In its default implementation it calls internally the function
+ * project to manifold. User classes can get away by simply
+ * implementing that method.
+ */
+ virtual
+ Point<spacedim>
+ get_new_point(const Quadrature<spacedim> &quad) const;
+
+ /**
+ * Given a point which lies close to the given manifold, it modifies
+ * it and projects it to manifold itself.
+ *
+ * This class is used by the default implementation of the function
+ * get_new_point(). It should be made pure virtual, but for
+ * historical reason, derived classes like Boundary<dim, spacedim>
+ * do not implement it. The default behavior of this class, however,
+ * is to throw an exception when called.
+ *
+ * If your manifold is simple, you could implement this function
+ * only, and the default behavior should work out of the box.
+ */
+ virtual
+ Point<spacedim> project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
+ const Point<spacedim> &candidate) const;
+
+ /**
+ * Backward compatibility interface. 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). The default implementation of this function
+ * passes its argument to the Manifolds::get_default_quadrature()
+ * function, and then calls the
+ * Manifold<dim,spacedim>::get_new_point() function. User derived
+ * classes can overload Manifold<dim,spacedim>::get_new_point() or
+ * Manifold<dim,spacedim>::project_to_surface(), which is called by
+ * the default implementation of
+ * Manifold<dim,spacedim>::get_new_point().
+ */
+ virtual
+ Point<spacedim>
+ get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
+
+ /**
+ * Backward compatibility interface. 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>.
+ *
+ * The default implementation of this function passes its argument
+ * to the Manifolds::get_default_quadrature() function, and then
+ * calls the Manifold<dim,spacedim>::get_new_point() function. User
+ * derived classes can overload
+ * Manifold<dim,spacedim>::get_new_point() or
+ * Manifold<dim,spacedim>::project_to_surface(), which is called by
+ * the default implementation of
+ * Manifold<dim,spacedim>::get_new_point().
+ */
+ virtual
+ Point<spacedim>
+ get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
+
+ /**
+ * Backward compatibility interface. Return the point which shall
+ * become the common point of the eight children of a hex in three
+ * or 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 all the bounding objects of the
+ * given @p hex are refined, so you may want to use the information
+ * provided by <tt>hex->quad(i)->line(j)->child(k)</tt>,
+ * <tt>i=0...5</tt>, <tt>j=0...3</tt>, <tt>k=0,1</tt>.
+ *
+ * The default implementation of this function passes its argument
+ * to the Manifolds::get_default_quadrature() function, and then
+ * calls the Manifold<dim,spacedim>::get_new_point() function. User
+ * derived classes can overload
+ * Manifold<dim,spacedim>::get_new_point() or
+ * Manifold<dim,spacedim>::project_to_surface(), which is called by
+ * the default implementation of
+ * Manifold<dim,spacedim>::get_new_point().
+ */
+ virtual
+ Point<spacedim>
+ get_new_point_on_hex (const typename Triangulation<dim,spacedim>::hex_iterator &hex) const;
+
+
+ /**
+ * Backward compatibility interface. Depending on <tt>dim=2</tt> or
+ * <tt>dim=3</tt> this function calls the get_new_point_on_line or
+ * the get_new_point_on_quad function. It throws an exception for
+ * <tt>dim=1</tt>. This wrapper allows dimension independent
+ * programming.
+ */
+ Point<spacedim>
+ get_new_point_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face) const;
+
+
+ /**
+ * Backward compatibility interface. Depending on <tt>dim=1</tt>,
+ * <tt>dim=2</tt> or <tt>dim=3</tt> this function calls the
+ * get_new_point_on_line, get_new_point_on_quad or the
+ * get_new_point_on_hex function. This wrapper allows dimension
+ * independent programming.
+ */
+ Point<spacedim>
+ get_new_point_on_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
+};
+
+
+/**
+ * Specialization of Manifold<dim,spacedim>, which represent a
+ * possibly periodic Euclidean space of dimension @p dim embedded in
+ * the Euclidean space of @p spacedim dimensions. The main
+ * characteristic of this Manifold is the fact that the function
+ * FlatManifold<dim,spacedim>::project_to_manifold() is the identity
+ * function.
+ *
+ * @ingroup manifold
+ *
+ * @author Luca Heltai, 2014
+ */
+template <int dim, int spacedim=dim>
+class FlatManifold: public Manifold<dim, spacedim>
+{
+public:
+ /**
+ * Default constructor.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). Special cases
+ * are taken into account, to ensure that the behavior is always as
+ * expected.
+ *
+ * Periodicity will be intended in the following way: the domain is
+ * considered to be the box contained in [Point<spacedim>(),
+ * periodicity) where the right extreme is excluded. If any of the
+ * components of this box has zero length, then no periodicity is
+ * computed in that direction. Whenever a function that tries to
+ * compute averages is called, an exception will be thrown if one of
+ * the points which you are using for the average lies outside the
+ * periodicity box. The return points are garanteed to lie in the
+ * perodicity box.
+ */
+ 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 calls 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 Quadrature<spacedim> &quad) const;
+
+
+ /**
+ * Project to FlatManifold. This is the identity function for flat,
+ * Euclidean spaces. 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 chartdim, which is part of a
+ * Manifold<dim,spacedim>. This object specializes a Manifold of
+ * dimension chartdim embedded in a manifold of dimension spacedim,
+ * for which you have explicit pull_back and push_forward
+ * transformations. This object only makes sense when chartdim <=
+ * dim, and the constructor throws an exception if this is not the
+ * case.
+ *
+ * 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 represents your manifold, i.e.,
+ * when your manifold \f$\mathcal{M}\f$ can be represented by a map
+ * \f[
+ * F: \mathcal{B} \subset R^{\text{chartdim}} \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{chartdim}}
+ * \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
+ * chartdim 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=dim, int chartdim=dim>
+class ManifoldChart: public Manifold<dim,spacedim>
+{
+public:
+ /**
+ * Constructor. The optional argument can be used to specify the
+ * periodicity of the chartdim-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<chartdim> periodicity=Point<chartdim>());
+
+ /**
+ * 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 Quadrature<spacedim> &quad) 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<chartdim>
+ 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<chartdim> &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,spacedim> sub_manifold;
+};
+
+
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
+#ifndef DOXYGEN
+
+template <>
+Point<1>
+Manifold<1,1>::
+get_new_point_on_face (const Triangulation<1,1>::face_iterator &) const;
+
+template <>
+Point<2>
+Manifold<1,2>::
+get_new_point_on_face (const Triangulation<1,2>::face_iterator &) const;
+
+
+template <>
+Point<3>
+Manifold<1,3>::
+get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const;
+
+
+template <>
+Point<1>
+Manifold<1,1>::
+get_new_point_on_quad (const Triangulation<1,1>::quad_iterator &) const;
+
+template <>
+Point<2>
+Manifold<1,2>::
+get_new_point_on_quad (const Triangulation<1,2>::quad_iterator &) const;
+
+
+template <>
+Point<3>
+Manifold<1,3>::
+get_new_point_on_quad (const Triangulation<1,3>::quad_iterator &) const;
+
+
+template <>
+Point<3>
+Manifold<3,3>::
+get_new_point_on_hex (const Triangulation<3,3>::hex_iterator &) const;
+
+#endif // DOXYGEN
+
+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
* around a given center point.
*
* @ingroup boundary
- * @author Wolfgang Bangerth, 1999, 2001, 2009, Ralf Hartmann, 2001, 2008
+ * @author Wolfgang Bangerth, 1999, 2001, 2009, Ralf Hartmann, 2001, 2008, Luca Heltai, 2014
*/
template <int dim, int spacedim=dim>
-class Boundary : public Subscriptor
+class Boundary : public FlatManifold<dim, spacedim>
{
public:
*/
virtual ~Boundary ();
- /**
- * 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<spacedim>
- get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const = 0;
-
- /**
- * 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<spacedim>
- get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
-
- /**
- * Depending on <tt>dim=2</tt> or <tt>dim=3</tt> this function calls the
- * get_new_point_on_line or the get_new_point_on_quad function. It throws an
- * exception for <tt>dim=1</tt>. This wrapper allows dimension independent
- * programming.
- */
- Point<spacedim>
- get_new_point_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face) const;
/**
* Return intermediate points on a line spaced according to the interior
#ifndef DOXYGEN
-template <>
-Point<1>
-Boundary<1,1>::
-get_new_point_on_face (const Triangulation<1,1>::face_iterator &) const;
-
template <>
void
Boundary<1,1>::
get_intermediate_points_on_face (const Triangulation<1,1>::face_iterator &,
std::vector<Point<1> > &) const;
-template <>
-Point<2>
-Boundary<1,2>::
-get_new_point_on_face (const Triangulation<1,2>::face_iterator &) const;
-
template <>
void
Boundary<1,2>::
get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &,
std::vector<Point<2> > &) const;
-
-
-template <>
-Point<3>
-Boundary<1,3>::
-get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const;
-
template <>
void
Boundary<1,3>::
grid_reordering.cc
grid_tools.cc
intergrid_map.cc
+ manifold.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
tria_accessor.inst.in
tria_boundary.inst.in
tria_boundary_lib.inst.in
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// 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/tensor.h>
+#include <deal.II/grid/manifold.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/fe/fe_q.h>
+#include <cmath>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace Manifolds {
+
+ void get_default_quadrature(const TriaIterator<CellAccessor<3, 3> >& obj,
+ Quadrature<3> &quad)
+ {
+ std::vector<Point<3> > sp;
+ std::vector<double> wp;
+
+ const int dim = 3;
+
+ int np = GeometryInfo<dim>::vertices_per_cell+
+ GeometryInfo<dim>::lines_per_cell+
+ GeometryInfo<dim>::faces_per_cell;
+ sp.resize(np);
+ wp.resize(np);
+ unsigned int j=0;
+ for(unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i, ++j)
+ {
+ sp[j] = obj->vertex(i);
+ wp[j] = 1.0/128.0;
+ }
+ for(unsigned int i=0; i<GeometryInfo<dim>::lines_per_cell; ++i, ++j)
+ {
+ sp[j] = (obj->line(i)->has_children() ? obj->line(i)->child(0)->vertex(1) :
+ obj->line(i)->center());
+ wp[j] = 7.0/192.0;
+ }
+ for(unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i, ++j)
+ {
+ sp[j] = (obj->face(i)->has_children() ? obj->face(i)->isotropic_child(0)->vertex(3) :
+ obj->face(i)->center());
+ wp[j] = 1.0/12.0;
+ }
+ quad = Quadrature<3>(sp,wp);
+ }
+
+
+ template <typename OBJECT, int spacedim>
+ void get_default_quadrature(const OBJECT& obj,
+ Quadrature<spacedim> &quad, bool with_laplace = false)
+ {
+ std::vector<Point<spacedim> > sp;
+ std::vector<double> wp;
+
+ 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[1] = .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;
+ }
+ quad = Quadrature<spacedim>(sp,wp);
+ }
+}
+
+using namespace Manifolds;
+
+/* -------------------------- Manifold --------------------- */
+
+
+template <int dim, int spacedim>
+Manifold<dim, spacedim>::~Manifold ()
+{}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+Manifold<dim, spacedim>::
+project_to_manifold (const std::vector<Point<spacedim> > &,
+ const Point<spacedim> &) const
+{
+ Assert (false, ExcPureFunctionCalled());
+ return Point<spacedim>();
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+Manifold<dim, spacedim>::
+get_new_point (const Quadrature<spacedim> &) const
+{
+ Assert (false, ExcPureFunctionCalled());
+ return Point<spacedim>();
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+Manifold<dim, spacedim>::
+get_new_point_on_line (const typename Triangulation<dim, spacedim>::line_iterator &line) const
+{
+ Quadrature<spacedim> quadrature;
+ get_default_quadrature(line, quadrature);
+ return get_new_point(quadrature);
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+Manifold<dim, spacedim>::
+get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
+{
+ Quadrature<spacedim> quadrature;
+ get_default_quadrature(quad, quadrature);
+ return get_new_point(quadrature);
+}
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+Manifold<dim,spacedim>::
+get_new_point_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face) const
+{
+ Assert (dim>1, ExcImpossibleInDim(dim));
+
+ switch (dim)
+ {
+ case 2:
+ return get_new_point_on_line (face);
+ case 3:
+ return get_new_point_on_quad (face);
+ }
+
+ return Point<spacedim>();
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+Manifold<dim,spacedim>::
+get_new_point_on_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
+{
+ switch (dim)
+ {
+ case 1:
+ return get_new_point_on_line (cell);
+ case 2:
+ return get_new_point_on_quad (cell);
+ case 3:
+ return get_new_point_on_hex (cell);
+ }
+
+ return Point<spacedim>();
+}
+
+
+template <>
+Point<1>
+Manifold<1,1>::
+get_new_point_on_face (const Triangulation<1,1>::face_iterator &) const
+{
+ Assert (false, ExcImpossibleInDim(1));
+ return Point<1>();
+}
+
+
+template <>
+Point<2>
+Manifold<1,2>::
+get_new_point_on_face (const Triangulation<1,2>::face_iterator &) const
+{
+ Assert (false, ExcImpossibleInDim(1));
+ return Point<2>();
+}
+
+
+
+template <>
+Point<3>
+Manifold<1,3>::
+get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const
+{
+ Assert (false, ExcImpossibleInDim(1));
+ return Point<3>();
+}
+
+
+template <>
+Point<1>
+Manifold<1,1>::
+get_new_point_on_quad (const Triangulation<1,1>::quad_iterator &) const
+{
+ Assert (false, ExcImpossibleInDim(1));
+ return Point<1>();
+}
+
+template <>
+Point<2>
+Manifold<1,2>::
+get_new_point_on_quad (const Triangulation<1,2>::quad_iterator &) const
+{
+ Assert (false, ExcImpossibleInDim(1));
+ return Point<2>();
+}
+
+
+template <>
+Point<3>
+Manifold<1,3>::
+get_new_point_on_quad (const Triangulation<1,3>::quad_iterator &) const
+{
+ Assert (false, ExcImpossibleInDim(1));
+ return Point<3>();
+}
+
+template <int dim, int spacedim>
+Point<spacedim>
+Manifold<dim, spacedim>::
+get_new_point_on_hex (const typename Triangulation<dim, spacedim>::hex_iterator &hex) const
+{
+ Assert (false, ExcImpossibleInDim(dim));
+ return Point<spacedim>();
+}
+
+template <>
+Point<3>
+Manifold<3,3>::
+get_new_point_on_hex (const typename Triangulation<3, 3>::hex_iterator &hex) const{
+ Quadrature<3> quadrature;
+ get_default_quadrature(hex, quadrature);
+ return get_new_point(quadrature);
+}
+
+
+/* -------------------------- FlatManifold --------------------- */
+
+
+template <int dim, int spacedim>
+FlatManifold<dim,spacedim>::FlatManifold (const Point<spacedim> periodicity) :
+ periodicity(periodicity)
+{}
+
+template <int dim, int spacedim>
+Point<spacedim>
+FlatManifold<dim, spacedim>::
+get_new_point (const Quadrature<spacedim> &quad) const
+{
+ const std::vector<Point<spacedim> > &surrounding_points = quad.get_points();
+ const std::vector<double> &weights = quad.get_weights();
+
+#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;
+ Point<spacedim> minP = periodicity;
+ bool check_period = (periodicity.norm() != 0);
+ if(check_period)
+ for(unsigned int i=0; i<surrounding_points.size(); ++i)
+ for(unsigned int d=0; d<spacedim; ++d) {
+ minP[d] = std::min(minP[d], surrounding_points[i][d]);
+ if(periodicity[d] > 0)
+ Assert(surrounding_points[i][d] < periodicity[d],
+ ExcMessage("One of the points does not lye into the periodic box! Bailing out."));
+ }
+
+ for(unsigned int i=0; i<surrounding_points.size(); ++i) {
+ dp = Point<spacedim>();
+ if(check_period) {
+ for(unsigned int d=0; d<spacedim; ++d)
+ if(periodicity[d] > 0)
+ dp[d] = ( (surrounding_points[i][d]-minP[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] = (p[d] < 0 ? p[d] + periodicity[d] : p[d]);
+
+ return project_to_manifold(surrounding_points, p);
+}
+
+template <int dim, int spacedim>
+Point<spacedim>
+FlatManifold<dim, spacedim>::project_to_manifold (const std::vector<Point<spacedim> > & vertices,
+ const Point<spacedim> &candidate) const
+{
+ return candidate;
+}
+
+
+/* -------------------------- ManifoldChart --------------------- */
+
+template <int dim, int spacedim, int chartdim>
+ManifoldChart<dim,spacedim,chartdim>::~ManifoldChart ()
+{}
+
+template <int dim, int spacedim, int chartdim>
+ManifoldChart<dim,spacedim,chartdim>::ManifoldChart (const Point<chartdim> periodicity):
+ sub_manifold(periodicity)
+{}
+
+
+template <int dim, int spacedim, int chartdim>
+Point<spacedim>
+ManifoldChart<dim,spacedim,chartdim>::
+get_new_point (const Quadrature<spacedim> &quad) const
+{
+ const std::vector<Point<spacedim> > &surrounding_points = quad.get_points();
+ const std::vector<double> &weights = quad.get_weights();
+ std::vector<Point<chartdim> > chart_points(surrounding_points.size());
+
+ for(unsigned int i=0; i<surrounding_points.size(); ++i)
+ chart_points[i] = pull_back(surrounding_points[i]);
+
+ Point<chartdim> 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$
+//
+// 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_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
+ {
+#if deal_II_dimension <= deal_II_space_dimension
+ template class Manifold<deal_II_dimension, deal_II_space_dimension>;
+ template class FlatManifold<deal_II_dimension, deal_II_space_dimension>;
+#endif
+ }
+
+
+
// average of the
// 8 surrounding
// vertices
- Point<spacedim> 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;
+ // Point<spacedim> new_point;
+ Quadrature<spacedim> quadrature;
+ Manifolds::get_default_quadrature(cell,quadrature);
+
+ // triangulation.vertices[next_unused_vertex] = new_point;
+ triangulation.vertices[next_unused_vertex] =
+ cell->get_manifold().get_new_point(quadrature);
// if the user_flag is set, i.e. if the
// cell is at the boundary, use a
// 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<dim>
- ::opposite_face[boundary_face])
- ->child(0)->vertex(1));
- }
+ {
+ std::vector<Point<spacedim> > ps(2);
+ std::vector<double> ws(2, 0.5);
+ ps[0] = cell->face(boundary_face)
+ ->child(0)->vertex(1);
+ ps[1] = cell->face(GeometryInfo<dim>
+ ::opposite_face[boundary_face])
+ ->child(0)->vertex(1);
+ Quadrature<spacedim> qs(ps,ws);
+ triangulation.vertices[next_unused_vertex]
+ = cell->get_manifold().get_new_point(qs);
+ }
+ }
}
else
{
triangulation.vertices[next_unused_vertex]
= quad->get_manifold().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;
-
- triangulation.vertices_used[next_unused_vertex] = true;
-
+ {
+ // 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;
+ Quadrature<spacedim> qs;
+ Manifolds::get_default_quadrature(quad, qs, true);
+ triangulation.vertices[next_unused_vertex] =
+ quad->get_manifold().get_new_point (qs);
+ }
+ triangulation.vertices_used[next_unused_vertex] = true;
// now that we created
// the right point, make
// up the four lines
// vertex at the center of
// the quad is weighted (see
// above)
- triangulation.vertices[next_unused_vertex] = Point<dim>();
- // first add corners of hex
- for (unsigned int vertex=0;
- vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
- triangulation.vertices[next_unused_vertex] += hex->vertex(vertex) / 128;
- // now add center of lines
- for (unsigned int line=0;
- line<GeometryInfo<dim>::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<GeometryInfo<dim>::faces_per_cell; ++face)
- triangulation.vertices[next_unused_vertex] += hex->face(face)->isotropic_child(0)->vertex(3) *
- 1./12.;
+
+ // triangulation.vertices[next_unused_vertex] = Point<dim>();
+ // // first add corners of hex
+ // for (unsigned int vertex=0;
+ // vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
+ // triangulation.vertices[next_unused_vertex] += hex->vertex(vertex) / 128;
+ // // now add center of lines
+ // for (unsigned int line=0;
+ // line<GeometryInfo<dim>::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<GeometryInfo<dim>::faces_per_cell; ++face)
+ // triangulation.vertices[next_unused_vertex] += hex->face(face)->isotropic_child(0)->vertex(3) *
+ // 1./12.;
+
+ triangulation.vertices[next_unused_vertex] =
+ hex->get_manifold().get_new_point_on_hex(hex);
// set the data of the
// six lines. first
{}
-
-template <int dim, int spacedim>
-Point<spacedim>
-Boundary<dim, spacedim>::
-get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &) const
-{
- Assert (false, ExcPureFunctionCalled());
- return Point<spacedim>();
-}
-
-
-
template <int dim, int spacedim>
void
Boundary<dim, spacedim>::
}
-
-template <int dim, int spacedim>
-Point<spacedim>
-Boundary<dim,spacedim>::
-get_new_point_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face) const
-{
- Assert (dim>1, ExcImpossibleInDim(dim));
-
- switch (dim)
- {
- case 2:
- return get_new_point_on_line (face);
- case 3:
- return get_new_point_on_quad (face);
- }
-
- return Point<spacedim>();
-}
-
-
template <int dim, int spacedim>
void
Boundary<dim,spacedim>::
}
-
-template <>
-Point<1>
-Boundary<1,1>::
-get_new_point_on_face (const Triangulation<1,1>::face_iterator &) const
-{
- Assert (false, ExcImpossibleInDim(1));
- return Point<1>();
-}
-
-
template <>
void
Boundary<1,1>::
}
-
-template <>
-Point<2>
-Boundary<1,2>::
-get_new_point_on_face (const Triangulation<1,2>::face_iterator &) const
-{
- Assert (false, ExcImpossibleInDim(1));
- return Point<2>();
-}
-
-
template <>
void
Boundary<1,2>::
}
-
-template <>
-Point<3>
-Boundary<1,3>::
-get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const
-{
- Assert (false, ExcImpossibleInDim(1));
- return Point<3>();
-}
-
-
template <>
void
Boundary<1,3>::
template <int dim, int spacedim>
Point<spacedim>
StraightBoundary<dim, spacedim>::
-get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &) const
+get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
{
- Assert (false, ExcImpossibleInDim(dim));
- return Point<spacedim>();
+ return FlatManifold<dim,spacedim>::get_new_point_on_quad(quad);
}
--- /dev/null
+//---------------------------- manifold_id_01.cc ---------------------------
+// Copyright (C) 2011, 2013 by the mathLab team.
+//
+// This file is subject to LGPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- flat_manifold_01.cc ---------------------------
+
+
+// Test that the flat manifold does what it should
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+
+// all include files you need here
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_out.h>
+
+// Helper function
+template <int dim, int spacedim>
+void test(unsigned int ref=1)
+{
+ deallog << "Testing dim=" << dim
+ << ", spacedim="<< spacedim << std::endl;
+
+ Triangulation<dim,spacedim> tria;
+ GridGenerator::hyper_cube (tria);
+ tria.refine_global(1);
+
+ typename Triangulation<dim,spacedim>::active_cell_iterator
+ cell;
+
+ for(cell=tria.begin_active(); cell!=tria.end(); ++cell) {
+
+ // check that FlatManifold returns the middle of the cell.
+ deallog << "Cell: " << cell << std::endl;
+ if(cell->get_manifold().get_new_point_on_cell(cell).distance(cell->center()) > 1e-6) {
+ deallog << "Default manifold: " << cell->get_manifold().get_new_point_on_cell(cell) << std::endl;
+ deallog << "Center of cell : " << cell->center() << std::endl;
+ } else {
+ deallog << "OK!" << std::endl;
+ }
+ }
+}
+
+int main ()
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<1,1>();
+ test<1,2>();
+ test<2,2>();
+ test<2,3>();
+ test<3,3>();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::Testing dim=1, spacedim=1
+DEAL::Cell: 1.0
+DEAL::OK!
+DEAL::Cell: 1.1
+DEAL::OK!
+DEAL::Testing dim=1, spacedim=2
+DEAL::Cell: 1.0
+DEAL::OK!
+DEAL::Cell: 1.1
+DEAL::OK!
+DEAL::Testing dim=2, spacedim=2
+DEAL::Cell: 1.0
+DEAL::OK!
+DEAL::Cell: 1.1
+DEAL::OK!
+DEAL::Cell: 1.2
+DEAL::OK!
+DEAL::Cell: 1.3
+DEAL::OK!
+DEAL::Testing dim=2, spacedim=3
+DEAL::Cell: 1.0
+DEAL::OK!
+DEAL::Cell: 1.1
+DEAL::OK!
+DEAL::Cell: 1.2
+DEAL::OK!
+DEAL::Cell: 1.3
+DEAL::OK!
+DEAL::Testing dim=3, spacedim=3
+DEAL::Cell: 1.0
+DEAL::OK!
+DEAL::Cell: 1.1
+DEAL::OK!
+DEAL::Cell: 1.2
+DEAL::OK!
+DEAL::Cell: 1.3
+DEAL::OK!
+DEAL::Cell: 1.4
+DEAL::OK!
+DEAL::Cell: 1.5
+DEAL::OK!
+DEAL::Cell: 1.6
+DEAL::OK!
+DEAL::Cell: 1.7
+DEAL::OK!
--- /dev/null
+//---------------------------- manifold_id_01.cc ---------------------------
+// Copyright (C) 2011, 2013 by the mathLab team.
+//
+// This file is subject to LGPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- flat_manifold_02.cc ---------------------------
+
+
+// Test that the flat manifold does what it should. This time on faces.
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+
+// all include files you need here
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_out.h>
+
+// Helper function
+template <int dim, int spacedim>
+void test(unsigned int ref=1)
+{
+ deallog << "Testing dim=" << dim
+ << ", spacedim="<< spacedim << std::endl;
+
+ Triangulation<dim,spacedim> tria;
+ GridGenerator::hyper_cube (tria);
+ tria.refine_global(1);
+
+ typename Triangulation<dim,spacedim>::active_cell_iterator
+ cell;
+
+ for(cell=tria.begin_active(); cell!=tria.end(); ++cell) {
+
+ // check that FlatManifold returns the middle of the cell.
+ deallog << "Cell: " << cell << std::endl;
+ for(unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f) {
+ const typename Triangulation<dim,spacedim>::face_iterator &face = cell->face(f);
+ if(face->get_manifold().get_new_point_on_face(face).distance(face->center()) > 1e-6) {
+ deallog << "Face : " << face << std::endl;
+ deallog << "Default manifold: " << cell->get_manifold().get_new_point_on_cell(cell) << std::endl;
+ deallog << "Center of cell : " << cell->center() << std::endl;
+ } else {
+ deallog << "Face " << face << " is OK!" << std::endl;
+ }
+ }
+ }
+}
+
+int main ()
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<2,2>();
+ test<2,3>();
+ test<3,3>();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::Testing dim=2, spacedim=2
+DEAL::Cell: 1.0
+DEAL::Face 6 is OK!
+DEAL::Face 12 is OK!
+DEAL::Face 4 is OK!
+DEAL::Face 14 is OK!
+DEAL::Cell: 1.1
+DEAL::Face 12 is OK!
+DEAL::Face 8 is OK!
+DEAL::Face 5 is OK!
+DEAL::Face 15 is OK!
+DEAL::Cell: 1.2
+DEAL::Face 7 is OK!
+DEAL::Face 13 is OK!
+DEAL::Face 14 is OK!
+DEAL::Face 10 is OK!
+DEAL::Cell: 1.3
+DEAL::Face 13 is OK!
+DEAL::Face 9 is OK!
+DEAL::Face 15 is OK!
+DEAL::Face 11 is OK!
+DEAL::Testing dim=2, spacedim=3
+DEAL::Cell: 1.0
+DEAL::Face 6 is OK!
+DEAL::Face 12 is OK!
+DEAL::Face 4 is OK!
+DEAL::Face 14 is OK!
+DEAL::Cell: 1.1
+DEAL::Face 12 is OK!
+DEAL::Face 8 is OK!
+DEAL::Face 5 is OK!
+DEAL::Face 15 is OK!
+DEAL::Cell: 1.2
+DEAL::Face 7 is OK!
+DEAL::Face 13 is OK!
+DEAL::Face 14 is OK!
+DEAL::Face 10 is OK!
+DEAL::Cell: 1.3
+DEAL::Face 13 is OK!
+DEAL::Face 9 is OK!
+DEAL::Face 15 is OK!
+DEAL::Face 11 is OK!
+DEAL::Testing dim=3, spacedim=3
+DEAL::Cell: 1.0
+DEAL::Face 14 is OK!
+DEAL::Face 41 is OK!
+DEAL::Face 6 is OK!
+DEAL::Face 37 is OK!
+DEAL::Face 10 is OK!
+DEAL::Face 33 is OK!
+DEAL::Cell: 1.1
+DEAL::Face 41 is OK!
+DEAL::Face 18 is OK!
+DEAL::Face 8 is OK!
+DEAL::Face 35 is OK!
+DEAL::Face 11 is OK!
+DEAL::Face 32 is OK!
+DEAL::Cell: 1.2
+DEAL::Face 15 is OK!
+DEAL::Face 40 is OK!
+DEAL::Face 37 is OK!
+DEAL::Face 22 is OK!
+DEAL::Face 12 is OK!
+DEAL::Face 31 is OK!
+DEAL::Cell: 1.3
+DEAL::Face 40 is OK!
+DEAL::Face 19 is OK!
+DEAL::Face 35 is OK!
+DEAL::Face 24 is OK!
+DEAL::Face 13 is OK!
+DEAL::Face 30 is OK!
+DEAL::Cell: 1.4
+DEAL::Face 16 is OK!
+DEAL::Face 39 is OK!
+DEAL::Face 7 is OK!
+DEAL::Face 36 is OK!
+DEAL::Face 33 is OK!
+DEAL::Face 26 is OK!
+DEAL::Cell: 1.5
+DEAL::Face 39 is OK!
+DEAL::Face 20 is OK!
+DEAL::Face 9 is OK!
+DEAL::Face 34 is OK!
+DEAL::Face 32 is OK!
+DEAL::Face 27 is OK!
+DEAL::Cell: 1.6
+DEAL::Face 17 is OK!
+DEAL::Face 38 is OK!
+DEAL::Face 36 is OK!
+DEAL::Face 23 is OK!
+DEAL::Face 31 is OK!
+DEAL::Face 28 is OK!
+DEAL::Cell: 1.7
+DEAL::Face 38 is OK!
+DEAL::Face 21 is OK!
+DEAL::Face 34 is OK!
+DEAL::Face 25 is OK!
+DEAL::Face 30 is OK!
+DEAL::Face 29 is OK!
--- /dev/null
+//---------------------------- manifold_id_01.cc ---------------------------
+// Copyright (C) 2011, 2013 by the mathLab team.
+//
+// This file is subject to LGPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- flat_manifold_03.cc ---------------------------
+
+
+// Test periodicity of FlatManifold
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+
+// all include files you need here
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_out.h>
+
+// Helper function
+template <int dim, int spacedim>
+void test(unsigned int ref=1)
+{
+ deallog << "Testing dim=" << dim
+ << ", spacedim="<< spacedim << std::endl;
+
+ Point<spacedim> periodicity;
+ periodicity[0] = 5.0;
+
+ FlatManifold<dim,spacedim> manifold(periodicity);
+
+ Quadrature<spacedim> quad;
+ std::vector<std::vector<Point<spacedim> > >ps(8,std::vector<Point<spacedim> >(2));
+ Point<spacedim> middle;
+ std::vector<double > ws(2, 0.5);
+
+ // Case 1: both points are close to left boundary of periodicity
+ ps[0][0][0] = 1;
+ ps[0][1][0] = 2;
+ // Case 2: same, with different order
+ ps[1][0][0] = 2;
+ ps[1][1][0] = 1;
+ // Case 3: one is close to left, one to right
+ ps[2][0][0] = 1;
+ ps[2][1][0] = 4;
+ // Case 3: same, opposite order
+ ps[3][0][0] = 4;
+ ps[3][1][0] = 1;
+ // Case 3: both close to right
+ ps[4][0][0] = 3;
+ ps[4][1][0] = 4;
+ // Case 3: same, opposite order
+ ps[5][0][0] = 4;
+ ps[5][1][0] = 3;
+ // Case 3: both close to middle
+ ps[6][0][0] = 2;
+ ps[6][1][0] = 3;
+ // Case 3: same, opposite order
+ ps[7][0][0] = 3;
+ ps[7][1][0] = 2;
+
+ for(unsigned int i=0; i<ps.size(); ++i) {
+ quad = Quadrature<spacedim>(ps[i],ws);
+ middle = manifold.get_new_point(quad);
+ deallog << "P0: " << ps[i][0] << " , P1: " << ps[i][1] << " , Middle: " << middle << std::endl;
+ }
+
+}
+
+int main ()
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<1,1>();
+ test<1,2>();
+ test<2,2>();
+ test<2,3>();
+ test<3,3>();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::Testing dim=1, spacedim=1
+DEAL::P0: 1.00000 , P1: 2.00000 , Middle: 1.50000
+DEAL::P0: 2.00000 , P1: 1.00000 , Middle: 1.50000
+DEAL::P0: 1.00000 , P1: 4.00000 , Middle: 0.00000
+DEAL::P0: 4.00000 , P1: 1.00000 , Middle: 0.00000
+DEAL::P0: 3.00000 , P1: 4.00000 , Middle: 3.50000
+DEAL::P0: 4.00000 , P1: 3.00000 , Middle: 3.50000
+DEAL::P0: 2.00000 , P1: 3.00000 , Middle: 2.50000
+DEAL::P0: 3.00000 , P1: 2.00000 , Middle: 2.50000
+DEAL::Testing dim=1, spacedim=2
+DEAL::P0: 1.00000 0.00000 , P1: 2.00000 0.00000 , Middle: 1.50000 0.00000
+DEAL::P0: 2.00000 0.00000 , P1: 1.00000 0.00000 , Middle: 1.50000 0.00000
+DEAL::P0: 1.00000 0.00000 , P1: 4.00000 0.00000 , Middle: 0.00000 0.00000
+DEAL::P0: 4.00000 0.00000 , P1: 1.00000 0.00000 , Middle: 0.00000 0.00000
+DEAL::P0: 3.00000 0.00000 , P1: 4.00000 0.00000 , Middle: 3.50000 0.00000
+DEAL::P0: 4.00000 0.00000 , P1: 3.00000 0.00000 , Middle: 3.50000 0.00000
+DEAL::P0: 2.00000 0.00000 , P1: 3.00000 0.00000 , Middle: 2.50000 0.00000
+DEAL::P0: 3.00000 0.00000 , P1: 2.00000 0.00000 , Middle: 2.50000 0.00000
+DEAL::Testing dim=2, spacedim=2
+DEAL::P0: 1.00000 0.00000 , P1: 2.00000 0.00000 , Middle: 1.50000 0.00000
+DEAL::P0: 2.00000 0.00000 , P1: 1.00000 0.00000 , Middle: 1.50000 0.00000
+DEAL::P0: 1.00000 0.00000 , P1: 4.00000 0.00000 , Middle: 0.00000 0.00000
+DEAL::P0: 4.00000 0.00000 , P1: 1.00000 0.00000 , Middle: 0.00000 0.00000
+DEAL::P0: 3.00000 0.00000 , P1: 4.00000 0.00000 , Middle: 3.50000 0.00000
+DEAL::P0: 4.00000 0.00000 , P1: 3.00000 0.00000 , Middle: 3.50000 0.00000
+DEAL::P0: 2.00000 0.00000 , P1: 3.00000 0.00000 , Middle: 2.50000 0.00000
+DEAL::P0: 3.00000 0.00000 , P1: 2.00000 0.00000 , Middle: 2.50000 0.00000
+DEAL::Testing dim=2, spacedim=3
+DEAL::P0: 1.00000 0.00000 0.00000 , P1: 2.00000 0.00000 0.00000 , Middle: 1.50000 0.00000 0.00000
+DEAL::P0: 2.00000 0.00000 0.00000 , P1: 1.00000 0.00000 0.00000 , Middle: 1.50000 0.00000 0.00000
+DEAL::P0: 1.00000 0.00000 0.00000 , P1: 4.00000 0.00000 0.00000 , Middle: 0.00000 0.00000 0.00000
+DEAL::P0: 4.00000 0.00000 0.00000 , P1: 1.00000 0.00000 0.00000 , Middle: 0.00000 0.00000 0.00000
+DEAL::P0: 3.00000 0.00000 0.00000 , P1: 4.00000 0.00000 0.00000 , Middle: 3.50000 0.00000 0.00000
+DEAL::P0: 4.00000 0.00000 0.00000 , P1: 3.00000 0.00000 0.00000 , Middle: 3.50000 0.00000 0.00000
+DEAL::P0: 2.00000 0.00000 0.00000 , P1: 3.00000 0.00000 0.00000 , Middle: 2.50000 0.00000 0.00000
+DEAL::P0: 3.00000 0.00000 0.00000 , P1: 2.00000 0.00000 0.00000 , Middle: 2.50000 0.00000 0.00000
+DEAL::Testing dim=3, spacedim=3
+DEAL::P0: 1.00000 0.00000 0.00000 , P1: 2.00000 0.00000 0.00000 , Middle: 1.50000 0.00000 0.00000
+DEAL::P0: 2.00000 0.00000 0.00000 , P1: 1.00000 0.00000 0.00000 , Middle: 1.50000 0.00000 0.00000
+DEAL::P0: 1.00000 0.00000 0.00000 , P1: 4.00000 0.00000 0.00000 , Middle: 0.00000 0.00000 0.00000
+DEAL::P0: 4.00000 0.00000 0.00000 , P1: 1.00000 0.00000 0.00000 , Middle: 0.00000 0.00000 0.00000
+DEAL::P0: 3.00000 0.00000 0.00000 , P1: 4.00000 0.00000 0.00000 , Middle: 3.50000 0.00000 0.00000
+DEAL::P0: 4.00000 0.00000 0.00000 , P1: 3.00000 0.00000 0.00000 , Middle: 3.50000 0.00000 0.00000
+DEAL::P0: 2.00000 0.00000 0.00000 , P1: 3.00000 0.00000 0.00000 , Middle: 2.50000 0.00000 0.00000
+DEAL::P0: 3.00000 0.00000 0.00000 , P1: 2.00000 0.00000 0.00000 , Middle: 2.50000 0.00000 0.00000
GridOut gridout;
gridout.write_msh(tria, deallog.get_file_stream());
+
}
int main ()
projects / dealii.git / commitdiff