From 3c0ff06ff22329154afb10e1aa8ba641fa47f177 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Fri, 25 Jul 2014 23:05:27 +0200 Subject: [PATCH] Added manifold, copying over the old concept, but retaining dim. Made Boundary derived from Manifold. 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, whose default has been moved to Manifold Fixed periodicity of FlatManifold, and added a test for it. Fixed comment on periodicity. Removed commented instantiations. --- include/deal.II/grid/manifold.h | 471 +++++++++++++++++++++++++ include/deal.II/grid/tria_boundary.h | 58 +-- source/grid/CMakeLists.txt | 2 + source/grid/manifold.cc | 395 +++++++++++++++++++++ source/grid/manifold.inst.in | 28 ++ source/grid/tria.cc | 165 ++++----- source/grid/tria_boundary.cc | 70 +--- tests/manifold/flat_manifold_01.cc | 69 ++++ tests/manifold/flat_manifold_01.output | 46 +++ tests/manifold/flat_manifold_02.cc | 71 ++++ tests/manifold/flat_manifold_02.output | 100 ++++++ tests/manifold/flat_manifold_03.cc | 92 +++++ tests/manifold/flat_manifold_03.output | 46 +++ tests/manifold/manifold_id_06.cc | 1 + 14 files changed, 1403 insertions(+), 211 deletions(-) create mode 100644 include/deal.II/grid/manifold.h create mode 100644 source/grid/manifold.cc create mode 100644 source/grid/manifold.inst.in create mode 100644 tests/manifold/flat_manifold_01.cc create mode 100644 tests/manifold/flat_manifold_01.output create mode 100644 tests/manifold/flat_manifold_02.cc create mode 100644 tests/manifold/flat_manifold_02.output create mode 100644 tests/manifold/flat_manifold_03.cc create mode 100644 tests/manifold/flat_manifold_03.output diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h new file mode 100644 index 0000000000..4b6de55b67 --- /dev/null +++ b/include/deal.II/grid/manifold.h @@ -0,0 +1,471 @@ +// --------------------------------------------------------------------- +// $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 +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +template class Triangulation; + + +/** We collect here some helper functions used in the + Manifold 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 >& 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 + void get_default_quadrature(const OBJECT& obj, + Quadrature &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 new_vertex = manifold.get_new_point (quadrature); + * ... + * @endcode + * @p quadrature is a Quadrature 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, which are all derived from + * FlatManifold, allowing old user codes to keep using + * their boundary descriptors as Manifold 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 +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 + get_new_point(const Quadrature &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 + * 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 project_to_manifold (const std::vector > &surrounding_points, + const Point &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::get_new_point() function. User derived + * classes can overload Manifold::get_new_point() or + * Manifold::project_to_surface(), which is called by + * the default implementation of + * Manifold::get_new_point(). + */ + virtual + Point + get_new_point_on_line (const typename Triangulation::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 + * quad->line(i)->child(j), i=0...3, j=0,1. + * + * The default implementation of this function passes its argument + * to the Manifolds::get_default_quadrature() function, and then + * calls the Manifold::get_new_point() function. User + * derived classes can overload + * Manifold::get_new_point() or + * Manifold::project_to_surface(), which is called by + * the default implementation of + * Manifold::get_new_point(). + */ + virtual + Point + get_new_point_on_quad (const typename Triangulation::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 hex->quad(i)->line(j)->child(k), + * i=0...5, j=0...3, k=0,1. + * + * The default implementation of this function passes its argument + * to the Manifolds::get_default_quadrature() function, and then + * calls the Manifold::get_new_point() function. User + * derived classes can overload + * Manifold::get_new_point() or + * Manifold::project_to_surface(), which is called by + * the default implementation of + * Manifold::get_new_point(). + */ + virtual + Point + get_new_point_on_hex (const typename Triangulation::hex_iterator &hex) const; + + + /** + * Backward compatibility interface. Depending on dim=2 or + * dim=3 this function calls the get_new_point_on_line or + * the get_new_point_on_quad function. It throws an exception for + * dim=1. This wrapper allows dimension independent + * programming. + */ + Point + get_new_point_on_face (const typename Triangulation::face_iterator &face) const; + + + /** + * Backward compatibility interface. Depending on dim=1, + * dim=2 or dim=3 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 + get_new_point_on_cell (const typename Triangulation::cell_iterator &cell) const; +}; + + +/** + * Specialization of Manifold, 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::project_to_manifold() is the identity + * function. + * + * @ingroup manifold + * + * @author Luca Heltai, 2014 + */ +template +class FlatManifold: public Manifold +{ +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(), + * 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 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 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 + get_new_point(const Quadrature &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 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 chartdim, which is part of a + * Manifold. 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 +class ManifoldChart: public Manifold +{ +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 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 Quadrature &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 + 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 + +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 diff --git a/include/deal.II/grid/tria_boundary.h b/include/deal.II/grid/tria_boundary.h index ca948550fc..9ee79ac5e7 100644 --- a/include/deal.II/grid/tria_boundary.h +++ b/include/deal.II/grid/tria_boundary.h @@ -26,6 +26,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -79,10 +80,10 @@ template class Triangulation; * 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 -class Boundary : public Subscriptor +class Boundary : public FlatManifold { public: @@ -105,42 +106,6 @@ 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 - get_new_point_on_line (const typename Triangulation::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 - * 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 - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; - - /** - * Depending on dim=2 or dim=3 this function calls the - * get_new_point_on_line or the get_new_point_on_quad function. It throws an - * exception for dim=1. This wrapper allows dimension independent - * programming. - */ - Point - get_new_point_on_face (const typename Triangulation::face_iterator &face) const; /** * Return intermediate points on a line spaced according to the interior @@ -464,35 +429,18 @@ public: #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 > &) 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 > &) const; - - -template <> -Point<3> -Boundary<1,3>:: -get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const; - template <> void Boundary<1,3>:: diff --git a/source/grid/CMakeLists.txt b/source/grid/CMakeLists.txt index fcac409a1b..1b074c9523 100644 --- a/source/grid/CMakeLists.txt +++ b/source/grid/CMakeLists.txt @@ -24,6 +24,7 @@ SET(_src grid_reordering.cc grid_tools.cc intergrid_map.cc + manifold.cc persistent_tria.cc tria_accessor.cc tria_boundary.cc @@ -41,6 +42,7 @@ SET(_inst grid_refinement.inst.in grid_tools.inst.in intergrid_map.inst.in + manifold.inst.in tria_accessor.inst.in tria_boundary.inst.in tria_boundary_lib.inst.in diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc new file mode 100644 index 0000000000..0d9205e259 --- /dev/null +++ b/source/grid/manifold.cc @@ -0,0 +1,395 @@ +// --------------------------------------------------------------------- +// $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 +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace Manifolds { + + void get_default_quadrature(const TriaIterator >& obj, + Quadrature<3> &quad) + { + std::vector > sp; + std::vector wp; + + const int dim = 3; + + int np = GeometryInfo::vertices_per_cell+ + GeometryInfo::lines_per_cell+ + GeometryInfo::faces_per_cell; + sp.resize(np); + wp.resize(np); + unsigned int j=0; + for(unsigned int i=0; i::vertices_per_cell; ++i, ++j) + { + sp[j] = obj->vertex(i); + wp[j] = 1.0/128.0; + } + for(unsigned int i=0; i::lines_per_cell; ++i, ++j) + { + sp[j] = (obj->line(i)->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::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 + void get_default_quadrature(const OBJECT& obj, + Quadrature &quad, bool with_laplace = false) + { + std::vector > sp; + std::vector 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(sp,wp); + } +} + +using namespace Manifolds; + +/* -------------------------- Manifold --------------------- */ + + +template +Manifold::~Manifold () +{} + + + +template +Point +Manifold:: +project_to_manifold (const std::vector > &, + const Point &) const +{ + Assert (false, ExcPureFunctionCalled()); + return Point(); +} + + + +template +Point +Manifold:: +get_new_point (const Quadrature &) const +{ + Assert (false, ExcPureFunctionCalled()); + return Point(); +} + + + +template +Point +Manifold:: +get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + Quadrature quadrature; + get_default_quadrature(line, quadrature); + return get_new_point(quadrature); +} + + + +template +Point +Manifold:: +get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +{ + Quadrature quadrature; + get_default_quadrature(quad, quadrature); + return get_new_point(quadrature); +} + + +template +Point +Manifold:: +get_new_point_on_face (const typename Triangulation::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(); +} + + + +template +Point +Manifold:: +get_new_point_on_cell (const typename Triangulation::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(); +} + + +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 +Point +Manifold:: +get_new_point_on_hex (const typename Triangulation::hex_iterator &hex) const +{ + Assert (false, ExcImpossibleInDim(dim)); + return Point(); +} + +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 +FlatManifold::FlatManifold (const Point periodicity) : + periodicity(periodicity) +{} + +template +Point +FlatManifold:: +get_new_point (const Quadrature &quad) const +{ + const std::vector > &surrounding_points = quad.get_points(); + const std::vector &weights = quad.get_weights(); + +#ifdef DEBUG + double sum=0; + for(unsigned int i=0; i p; + Point dp; + Point minP = periodicity; + bool check_period = (periodicity.norm() != 0); + if(check_period) + for(unsigned int i=0; i 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(); + if(check_period) { + for(unsigned int d=0; 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 0) + p[d] = (p[d] < 0 ? p[d] + periodicity[d] : p[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 Quadrature &quad) const +{ + const std::vector > &surrounding_points = quad.get_points(); + const std::vector &weights = quad.get_weights(); + 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/source/grid/manifold.inst.in b/source/grid/manifold.inst.in new file mode 100644 index 0000000000..892f3f92bf --- /dev/null +++ b/source/grid/manifold.inst.in @@ -0,0 +1,28 @@ +// --------------------------------------------------------------------- +// $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; + template class FlatManifold; +#endif + } + + + diff --git a/source/grid/tria.cc b/source/grid/tria.cc index eeeef5adc9..630c347f79 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -3961,12 +3961,13 @@ namespace internal // average of the // 8 surrounding // vertices - Point new_point; - for (unsigned int i=0; i<8; ++i) - new_point += triangulation.vertices[new_vertices[i]]; - new_point /= 8.0; - - triangulation.vertices[next_unused_vertex] = new_point; + // Point new_point; + Quadrature 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 @@ -4011,13 +4012,19 @@ namespace internal // 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 - ::opposite_face[boundary_face]) - ->child(0)->vertex(1)); - } + { + std::vector > ps(2); + std::vector ws(2, 0.5); + ps[0] = cell->face(boundary_face) + ->child(0)->vertex(1); + ps[1] = cell->face(GeometryInfo + ::opposite_face[boundary_face]) + ->child(0)->vertex(1); + Quadrature qs(ps,ws); + triangulation.vertices[next_unused_vertex] + = cell->get_manifold().get_new_point(qs); + } + } } else { @@ -6096,61 +6103,41 @@ namespace internal 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 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 @@ -8413,26 +8400,28 @@ namespace internal // vertex at the center of // the quad is weighted (see // above) - triangulation.vertices[next_unused_vertex] = Point(); - // first add corners of hex - for (unsigned int vertex=0; - vertex::vertices_per_cell; ++vertex) - triangulation.vertices[next_unused_vertex] += hex->vertex(vertex) / 128; - // now add center of lines - for (unsigned int line=0; - line::lines_per_cell; ++line) - triangulation.vertices[next_unused_vertex] += hex->line(line)->child(0)->vertex(1) * - 7./192.; - // finally add centers of - // faces. note that vertex 3 - // of child 0 is an invariant - // with respect to the face - // orientation, flip and - // rotation - for (unsigned int face=0; - face::faces_per_cell; ++face) - triangulation.vertices[next_unused_vertex] += hex->face(face)->isotropic_child(0)->vertex(3) * - 1./12.; + + // triangulation.vertices[next_unused_vertex] = Point(); + // // first add corners of hex + // for (unsigned int vertex=0; + // vertex::vertices_per_cell; ++vertex) + // triangulation.vertices[next_unused_vertex] += hex->vertex(vertex) / 128; + // // now add center of lines + // for (unsigned int line=0; + // line::lines_per_cell; ++line) + // triangulation.vertices[next_unused_vertex] += hex->line(line)->child(0)->vertex(1) * + // 7./192.; + // // finally add centers of faces. note that + // // vertex 3 of child 0 is an invariant with + // // respect to the face orientation, flip and + // // rotation + // for (unsigned int face=0; + // face::faces_per_cell; ++face) + // triangulation.vertices[next_unused_vertex] += hex->face(face)->isotropic_child(0)->vertex(3) * + // 1./12.; + + triangulation.vertices[next_unused_vertex] = + hex->get_manifold().get_new_point_on_hex(hex); // set the data of the // six lines. first diff --git a/source/grid/tria_boundary.cc b/source/grid/tria_boundary.cc index 098f82b776..7024425f42 100644 --- a/source/grid/tria_boundary.cc +++ b/source/grid/tria_boundary.cc @@ -34,18 +34,6 @@ Boundary::~Boundary () {} - -template -Point -Boundary:: -get_new_point_on_quad (const typename Triangulation::quad_iterator &) const -{ - Assert (false, ExcPureFunctionCalled()); - return Point(); -} - - - template void Boundary:: @@ -67,26 +55,6 @@ get_intermediate_points_on_quad (const typename Triangulation::qu } - -template -Point -Boundary:: -get_new_point_on_face (const typename Triangulation::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(); -} - - template void Boundary:: @@ -109,17 +77,6 @@ get_intermediate_points_on_face (const typename Triangulation::fac } - -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>:: @@ -130,17 +87,6 @@ get_intermediate_points_on_face (const Triangulation<1,1>::face_iterator &, } - -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>:: @@ -151,17 +97,6 @@ get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &, } - -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>:: @@ -375,10 +310,9 @@ namespace template Point StraightBoundary:: -get_new_point_on_quad (const typename Triangulation::quad_iterator &) const +get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const { - Assert (false, ExcImpossibleInDim(dim)); - return Point(); + return FlatManifold::get_new_point_on_quad(quad); } diff --git a/tests/manifold/flat_manifold_01.cc b/tests/manifold/flat_manifold_01.cc new file mode 100644 index 0000000000..80dea95140 --- /dev/null +++ b/tests/manifold/flat_manifold_01.cc @@ -0,0 +1,69 @@ +//---------------------------- 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 +#include + + +// all include files you need here +#include +#include +#include +#include +#include +#include + +// Helper function +template +void test(unsigned int ref=1) +{ + deallog << "Testing dim=" << dim + << ", spacedim="<< spacedim << std::endl; + + Triangulation tria; + GridGenerator::hyper_cube (tria); + tria.refine_global(1); + + typename Triangulation::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; +} + diff --git a/tests/manifold/flat_manifold_01.output b/tests/manifold/flat_manifold_01.output new file mode 100644 index 0000000000..d86174d6fc --- /dev/null +++ b/tests/manifold/flat_manifold_01.output @@ -0,0 +1,46 @@ + +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! diff --git a/tests/manifold/flat_manifold_02.cc b/tests/manifold/flat_manifold_02.cc new file mode 100644 index 0000000000..6884df3eaf --- /dev/null +++ b/tests/manifold/flat_manifold_02.cc @@ -0,0 +1,71 @@ +//---------------------------- 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 +#include + + +// all include files you need here +#include +#include +#include +#include +#include +#include + +// Helper function +template +void test(unsigned int ref=1) +{ + deallog << "Testing dim=" << dim + << ", spacedim="<< spacedim << std::endl; + + Triangulation tria; + GridGenerator::hyper_cube (tria); + tria.refine_global(1); + + typename Triangulation::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::faces_per_cell; ++f) { + const typename Triangulation::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; +} + diff --git a/tests/manifold/flat_manifold_02.output b/tests/manifold/flat_manifold_02.output new file mode 100644 index 0000000000..faa2760645 --- /dev/null +++ b/tests/manifold/flat_manifold_02.output @@ -0,0 +1,100 @@ + +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! diff --git a/tests/manifold/flat_manifold_03.cc b/tests/manifold/flat_manifold_03.cc new file mode 100644 index 0000000000..636a668106 --- /dev/null +++ b/tests/manifold/flat_manifold_03.cc @@ -0,0 +1,92 @@ +//---------------------------- 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 +#include + + +// all include files you need here +#include +#include +#include +#include +#include +#include + +// Helper function +template +void test(unsigned int ref=1) +{ + deallog << "Testing dim=" << dim + << ", spacedim="<< spacedim << std::endl; + + Point periodicity; + periodicity[0] = 5.0; + + FlatManifold manifold(periodicity); + + Quadrature quad; + std::vector > >ps(8,std::vector >(2)); + Point middle; + std::vector 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[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; +} + diff --git a/tests/manifold/flat_manifold_03.output b/tests/manifold/flat_manifold_03.output new file mode 100644 index 0000000000..2d84eb5c6b --- /dev/null +++ b/tests/manifold/flat_manifold_03.output @@ -0,0 +1,46 @@ + +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 diff --git a/tests/manifold/manifold_id_06.cc b/tests/manifold/manifold_id_06.cc index 0a96f74a0a..9378c5f4fa 100644 --- a/tests/manifold/manifold_id_06.cc +++ b/tests/manifold/manifold_id_06.cc @@ -58,6 +58,7 @@ void test(unsigned int ref=1) GridOut gridout; gridout.write_msh(tria, deallog.get_file_stream()); + } int main () -- 2.39.5