From df03f960fd83892f0b2cda70b3d3a167cc7d2033 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Mon, 4 Aug 2014 15:26:23 +0200 Subject: [PATCH] Fixed spelling and indentation. Replied to most comments by Wolfgang. Fixed constness of sub_manifold in ManifoldChart, fixed a few comments. Fixed constness in a few places. Removed references for simple doubles, and removed empty lines. Added back references to quad when they are arguments of get_new_point Made get_default_quadrature query itself the manifold when vertices are not available on subobjects. Fixed instatiations for invalid Triangulations. Fixed last reference, and failing test. Some changes to the documentation. --- doc/doxygen/headers/boundary.h | 26 ++- doc/doxygen/headers/manifold.h | 35 +-- doc/news/changes.h | 34 +++ include/deal.II/grid/manifold.h | 167 ++++++++++----- include/deal.II/grid/manifold_lib.h | 4 +- include/deal.II/grid/tria.h | 9 +- source/grid/manifold.cc | 101 +++------ source/grid/manifold_lib.cc | 12 +- source/grid/tria.cc | 37 +++- tests/manifold/manifold_id_06.cc | 1 - tests/manifold/spherical_manifold_01.cc | 2 +- tests/manifold/spherical_manifold_01.output | 226 ++++++++++---------- 12 files changed, 371 insertions(+), 283 deletions(-) diff --git a/doc/doxygen/headers/boundary.h b/doc/doxygen/headers/boundary.h index facdcfc1c3..ab51f0f49e 100644 --- a/doc/doxygen/headers/boundary.h +++ b/doc/doxygen/headers/boundary.h @@ -37,12 +37,12 @@ * derived from the Manifold and Boundary base classes describe the * geometry of a domain. One can then attach an object of a class * derived from this base classes to the Triangulation object using - * the Triangulation::set_boundary() (the old, deprecated function) - * or Triangulation::set_manifold() functions, and the Triangulation + * the Triangulation::set_boundary() or + * Triangulation::set_manifold() functions, and the Triangulation * will ask the manifold object where a new vertex should be located * upon mesh refinement. Several classes already exist to support * the boundary of the most common geometries, e.g., - * CylinderBoundary, HyperBallBoundary, or HyperShellBoundary. + * CylinderBoundary, HyperBallBoundary, or HyperShellBoundary. * *
  • Integration: When using higher order finite element methods, * it is often necessary to compute cell terms (like cell @@ -96,15 +96,19 @@ * The behavior of the Triangulation class w.r.t. geometry * descriptions is the following: Triangulation::set_boundary() and * Triangulation::set_manifold() do the exact same thing: they attach - * a manifold descriptor to the specified id. The first function is - * maintained for backward compatibility, but will be deprecated in - * future releases. + * a manifold descriptor to the specified id. The first function + * expects a Boundary descriptor (which is a specialization of a + * Manifold description) and is provided mainly for backward + * compatibility, while the second class expects a Manifold + * descriptor. Notice that the Triangulation class only uses the + * Manifold interface, and you could describe both the interior and + * the boundary of the domain using the same object. The additional + * information contained in the Boundary interface is related to the + * computation of the exact normals. * - * Whenever a new vertex is needed in an object that supports - * geometry dependent refinement (i.e., a face on the boundary, an - * interior face, or a cell in codimension one), the Triangulation - * queries the manifold_id of the object which needs refinement. If the - * manifold_id is set to numbers::invalid_manifold_id, then the + * Whenever a new vertex is needed in an object, the Triangulation + * queries the manifold_id of the object which needs refinement. If + * the manifold_id is set to numbers::invalid_manifold_id, then the * Triangulation queries the boundary_id (if the face is on the * boundary) or the material_id (if the Triangulation is of * codimension one and the object is a cell). If the previous queries diff --git a/doc/doxygen/headers/manifold.h b/doc/doxygen/headers/manifold.h index 7752907a6e..cd7ac68106 100644 --- a/doc/doxygen/headers/manifold.h +++ b/doc/doxygen/headers/manifold.h @@ -60,21 +60,32 @@ * * * - * In deal.II, a Manifold is seen as a collection of points, together with a - * notion of distance between points (on the manifold). New points are - * typically obtained by providing a local coordinate system on the manifold, - * identifying existing points in the local coordinate system (pulling them - * back using the local map to obtain their local coordinates), find the new - * point in the local coordinate system by weighted sums of the existing - * points, and transforming back the point in the real space (pushing it - * forward using the local map). (However, this pull back and push forward - * only happens in classes that implement manifolds, and your own class can - * also choose other ways of figuring out where a new point should be - * located.) + * In deal.II, a Manifold is seen as a collection of points, together + * with a notion of distance between points (on the manifold). New + * points are typically obtained by providing a local coordinate + * system on the manifold, identifying existing points in the local + * coordinate system (pulling them back using the local map to obtain + * their local coordinates), find the new point in the local + * coordinate system by weighted sums of the existing points, and + * transforming back the point in the real space (pushing it forward + * using the local map). The main class that implements this mechanism + * is the ManifoldChart class, and this is the class that users will + * likely overload for complex geometries. * * While this process is non trivial in most cases of interest, for * most of the trivial geometries, like cylinders, spheres or shells, - * we provide reasonable implementations. + * we provide reasonable implementations. + * + * The Boundary of a Triangulation is a special case of Manifold, for + * which additional informations can be useful in user codes, such as + * normal to surfaces or to curves. If your coarse mesh is reasonably + * shaped, you might be interested in only attaching a manifold + * description to boundary portion of your domain. This can be done + * using the set_boundary class, which take as arguments a Boundary + * object (derived from Manifold). Notice that Triangulation uses only + * the Manifold interface, not the Boundary interface. Other tools, + * however, might need to compute exact normals at quadrature points, + * and therefore a wrapper to query Boundary objects is provided. * * @see @ref GlossManifoldIndicator "Glossary entry on manifold * indicators" diff --git a/doc/news/changes.h b/doc/news/changes.h index 272bce21d2..0e92c848cb 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -97,6 +97,40 @@ inconvenience this causes.
      + +
    1. New: Added support for curved interior cells for all Triangulation + dimensions. +
      + A new Manifold class was introduced which only contains the + interface needed by Triangulation to refine objects, leaving all boundary + related functions in the class Boundary, which was made + derived from Manifold.
      + This new construction allows for curved interior cells, and custom refinement + strategies. At the moment the following Manifolds are supported: +
        +
      • FlatManifold: This class replaces the old + StraightBoundary, and it adds support for periodic + manifolds. This is the simplest class one can use to create new Manifold classes; +
        • +
      • ManifoldChart: This is one of the most general Manifold + one can think of. The user can overload the functions ManifoldChart::pull_back() and + ManifoldChart::push_forward(), to obtain a very general curved geometry, following + concepts typical of elasticity; +
        • +
      • SphericalManifold: A simple implementation of spherical coordinates + transformations. This manifold allows hyper shells with curved interior cells which + follow the natural shape of the shell; +
        • +
      +
      + The functions + Triangulation::set_boundary() and Triangulation::get_boundary() can still be used to + set and get Boundary objects instead of Manifold ones. For the get function, an exception + is thrown if a conversion to a valid Boundary class cannot be made on the fly. +
      + (Luca Heltai, 2014/08/06) +
      2. +
    2. Ported: The build system now supports CMake 3.0.
      (Matthias Maier, 2014/07/15) diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h index 176e10ff2a..3a8f51d855 100644 --- a/include/deal.II/grid/manifold.h +++ b/include/deal.II/grid/manifold.h @@ -32,26 +32,31 @@ DEAL_II_NAMESPACE_OPEN template class Triangulation; -/** We collect here some helper functions used in the - Manifold classes. +/** + * 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); +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. + */ + Quadrature<3> + get_default_quadrature(const TriaIterator >& hex); - /** 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.*/ + /** + * Given a general mesh iterator, construct a quadrature with the + * Laplace weights 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); + Quadrature + get_default_quadrature(const OBJECT& obj, bool with_laplace = false); } @@ -68,26 +73,17 @@ namespace Manifolds { * @endcode * @p quadrature is a Quadrature object, which contains a * collection of points in @p spacedim dimension, and a collection of - * weights. + * weights (Note that unlike almost all other cases in the library, + * we here interpret the points in the quadrature object to be in + * real space, not on the reference cell.) * - * Internally, the get_new_point() function calls the + * Internaly, the get_new_point() function calls the * project_to_manifold() function after computing the weighted - * average of the quadrature poitns. This allows end users to only + * average of the quadrature points. This allows end users to only * overload project_to_manifold() for simple situations. * * 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. - * - * The default behavior of these backward compatible interfaces is - * to construct a Quadrature object containting the - * vertices, midpoints of lines, and midpoints of quads with the - * correct weight, and call get_new_point() with this quadrature. If - * you need finer tuning for lines, quads or hexes, you can overload - * any of the get_new_point_on_* functions. + * overloaded. * * FlatManifold is the specialization from which StraigthBoundary is * derived, where the project_to_manifold() function is the identity. @@ -114,7 +110,7 @@ public: * points together with appropriate weights. * * In its default implementation it calls internally the function - * project to manifold. User classes can get away by simply + * project_to_manifold. User classes can get away by simply * implementing that method. */ virtual @@ -248,10 +244,10 @@ class FlatManifold: public Manifold { public: /** - * Default constructor.The optional argument can be used to specify + * 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 + * per direction). A periodicity value of zero means that along that + * direction there is no periodicity. By default no periodicity is * assumed. * * Periodicity affects the way a middle point is computed. It is @@ -288,10 +284,19 @@ public: * 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. + * trivial) perturbations of a flat manifold. This is the case + * whenever the coarse mesh is a decent approximation of the + * manifold geometry. In this case, the middle point of a cell is + * close to true middle point of the manifold, and a projection may + * suffice. + * + * 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. You will have good approximations also with large + * deformations, as long as in the coarsest mesh size you are trying + * to refine, the middle point is not too far from the manifold mid + * point, i.e., as long as the coarse mesh size is small enough. */ virtual Point get_new_point(const Quadrature &quad) const; @@ -343,12 +348,10 @@ private: * 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. + * transformations. * * This is an helper class which is useful when you have an explicit - * map from an Euclidean space of dimension dim to an Euclidean + * map from an Euclidean space of dimension chartdim 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[ @@ -374,6 +377,10 @@ private: * Derived classes are required to implement the push_forward() and * the pull_back() methods. * + * Notice that the dimenisions #chartdim and #spacedim can be + * arbitrary, as long as the transformation from $\mathcal{B}$ to + * $\mathcal{M}$ is invertible. + * * @ingroup manifold * * @author Luca Heltai, 2013 @@ -385,8 +392,8 @@ 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 + * direction). A periodicity value of zero means that along that + * direction there is no periodicity. By default no periodicity is * assumed. * * Periodicity affects the way a middle point is computed. It is @@ -415,7 +422,7 @@ public: get_new_point(const Quadrature &quad) const; /** - * Pull back the given point in spacedim to the Euclidean dim + * Pull back the given point in spacedim to the Euclidean chartdim * dimensional space. * * Refer to the general documentation of this class for more @@ -425,7 +432,7 @@ public: pull_back(const Point &space_point) const = 0; /** - * Given a point in the dim dimensianal Euclidean space, this + * Given a point in the chartdim dimensional Euclidean space, this * method returns a point on the manifold embedded in the spacedim * Euclidean space. * @@ -440,7 +447,7 @@ public: * The sub_manifold object is used to compute the average of the * points in the chart coordinates system. */ - FlatManifold sub_manifold; + const FlatManifold sub_manifold; }; @@ -489,6 +496,68 @@ Point<3> Manifold<3,3>:: get_new_point_on_hex (const Triangulation<3,3>::hex_iterator &) const; +/*---Templated functions---*/ + +namespace Manifolds { + + template + Quadrature + get_default_quadrature(const OBJECT& obj, bool with_laplace) + { + std::vector > sp; + std::vector wp; + + const int dim = OBJECT::AccessorType::structure_dimension; + + // 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 + 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)->get_manifold().get_new_point_on_line(obj->line(0)); + sp[5] = obj->line(1)->has_children() ? + obj->line(1)->child(0)->vertex(1) : + obj->line(1)->get_manifold().get_new_point_on_line(obj->line(1)); + sp[6] = obj->line(2)->has_children() ? + obj->line(2)->child(0)->vertex(1) : + obj->line(2)->get_manifold().get_new_point_on_line(obj->line(2)); + sp[7] = obj->line(3)->has_children() ? + obj->line(3)->child(0)->vertex(1) : + obj->line(3)->get_manifold().get_new_point_on_line(obj->line(3)); + 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; + } + return Quadrature(sp,wp); + } +} + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index d170e435f2..654c99ae05 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -40,7 +40,7 @@ DEAL_II_NAMESPACE_OPEN * The two dimensional implementation of this class works by * transforming points to spherical coordinates, taking the average in * that coordinate system, and then transforming back the point to - * cartesian coordinates. For the three dimensional case, we use a + * Cartesian coordinates. For the three dimensional case, we use a * simpler approach: we take the average of the norm of the points, * and use this value to shift the average point along the radial * direction. In order for this manifold to work correctly, it cannot @@ -59,7 +59,7 @@ public: /** * The Constructor takes the center of the spherical coordinates * system. This class uses the pull_back and push_forward mechanism - * to transform from cartesian to spherical coordinate systems, + * to transform from Cartesian to spherical coordinate systems, * taking into account the periodicity of base Manifold in two * dimensions, while in three dimensions it takes the middle point, * and project it along the radius using the average radius of the diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h index f2dabfed9e..df114567ea 100644 --- a/include/deal.II/grid/tria.h +++ b/include/deal.II/grid/tria.h @@ -3647,8 +3647,13 @@ template <> unsigned int Triangulation<1,3>::max_adjacent_cells () const; // ------------------------------------------------------------------- - - +// Explicit invalid things... +template <> +const Manifold<2,1> & Triangulation<2, 1>::get_manifold(const types::manifold_id) const; +template <> +const Manifold<3,1> & Triangulation<3, 1>::get_manifold(const types::manifold_id) const; +template <> +const Manifold<3,2> & Triangulation<3, 2>::get_manifold(const types::manifold_id) const; #endif // DOXYGEN diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 0575eb80d0..6707cb27ef 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -24,23 +24,29 @@ DEAL_II_NAMESPACE_OPEN - namespace Manifolds { - void get_default_quadrature(const TriaIterator >& obj, - Quadrature<3> &quad) + Quadrature<3> + get_default_quadrature(const TriaIterator >& obj) { std::vector > sp; std::vector wp; const int dim = 3; - int np = GeometryInfo::vertices_per_cell+ - GeometryInfo::lines_per_cell+ - GeometryInfo::faces_per_cell; + const unsigned int np = + GeometryInfo::vertices_per_cell+ + GeometryInfo::lines_per_cell+ + GeometryInfo::faces_per_cell; + sp.resize(np); wp.resize(np); unsigned int j=0; + + // note that the exact weights are chosen such as to minimize the + // distortion of the eight new hexes from the optimal shape; their + // derivation and values is copied over from the + // @p{MappingQ::set_laplace_on_vector} function for(unsigned int i=0; i::vertices_per_cell; ++i, ++j) { sp[j] = obj->vertex(i); @@ -49,71 +55,18 @@ namespace Manifolds { 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()); + obj->line(i)->get_manifold().get_new_point_on_line(obj->line(i))); 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()); + obj->face(i)->get_manifold().get_new_point_on_face(obj->face(i))); wp[j] = 1.0/12.0; } - quad = Quadrature<3>(sp,wp); - } - - - template - void get_default_quadrature(const OBJECT& obj, - Quadrature &quad, bool with_laplace) - { - 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); + return Quadrature<3>(sp,wp); } + } using namespace Manifolds; @@ -167,9 +120,9 @@ 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); + return get_new_point + (get_default_quadrature::line_iterator, + spacedim>(line, false)); } @@ -179,9 +132,9 @@ 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); + return get_new_point + (get_default_quadrature::quad_iterator, + spacedim>(quad, false)); } @@ -296,9 +249,7 @@ 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); + return get_new_point(get_default_quadrature(hex)); } @@ -331,7 +282,7 @@ get_new_point (const Quadrature &quad) const Point p; Point dp; Point minP = periodicity; - bool check_period = (periodicity.norm() > tolerance); + const bool check_period = (periodicity.norm() > tolerance); if(check_period) for(unsigned int i=0; i &quad) const for(unsigned int i=0; i chart_quad(chart_points, weights); - Point p_chart = sub_manifold.get_new_point(chart_quad); + const Quadrature chart_quad(chart_points, weights); + const Point p_chart = sub_manifold.get_new_point(chart_quad); return push_forward(p_chart); } diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index e041d06899..5f0483eee6 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -70,8 +70,8 @@ Point SphericalManifold::push_forward(const Point &spherical_point) const { Assert(spherical_point[0] >=0.0, ExcMessage("Negative radius for given point.")); - const double &rho = spherical_point[0]; - const double &theta = spherical_point[1]; + const double rho = spherical_point[0]; + const double theta = spherical_point[1]; Point p; if(rho > 1e-10) @@ -99,8 +99,8 @@ Point SphericalManifold::pull_back(const Point &space_point) const { const Point R = space_point-center; const double rho = R.norm(); - const double &x = R[0]; - const double &y = R[1]; + const double x = R[0]; + const double y = R[1]; Point p; p[0] = rho; @@ -113,13 +113,11 @@ SphericalManifold::pull_back(const Point &space_point) c break; case 3: { - const double &z = R[2]; + const double z = R[2]; p[2] = atan2(y,x); // phi if(p[2] < 0) p[2] += 2*numbers::PI; // phi is periodic p[1] = atan2(sqrt(x*x+y*y),z); // theta - // if(p[1] < 0) - // p[1] += 2*numbers::PI; } break; default: diff --git a/source/grid/tria.cc b/source/grid/tria.cc index afe95637c7..8169c5952a 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -3954,15 +3954,12 @@ namespace internal // boundary object if (dim == spacedim) { - // in a first step, compute the location of central - // vertex as the average of the 8 surrounding - // vertices, using Manifolds helper functions - 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); + cell->get_manifold().get_new_point + (Manifolds::get_default_quadrature + ::active_cell_iterator, + spacedim>(cell)); // if the user_flag is set, i.e. if the // cell is at the boundary, use a @@ -5828,10 +5825,11 @@ namespace internal // derivation and values is copied over from // the @p{MappingQ::set_laplace_on_vector} // function - Quadrature qs; - Manifolds::get_default_quadrature(quad, qs, true); triangulation.vertices[next_unused_vertex] = - quad->get_manifold().get_new_point (qs); + quad->get_manifold().get_new_point + (Manifolds::get_default_quadrature + ::quad_iterator, + spacedim>(quad, true)); } triangulation.vertices_used[next_unused_vertex] = true; // now that we created the right point, make up @@ -13291,7 +13289,26 @@ Triangulation::remove_refinement_listener (RefinementListener &li refinement_listener_map.erase (refinement_listener_map.find (&listener)); } +template <> +const Manifold<2,1> & Triangulation<2, 1>::get_manifold(const types::manifold_id) const { + Assert(false, ExcImpossibleInDim(1)); + static FlatManifold<2,1> flat; + return flat; +} + +template <> +const Manifold<3,1> & Triangulation<3, 1>::get_manifold(const types::manifold_id) const { + Assert(false, ExcImpossibleInDim(1)); + static FlatManifold<3,1> flat; + return flat; +} +template <> +const Manifold<3,2> & Triangulation<3, 2>::get_manifold(const types::manifold_id) const { + Assert(false, ExcImpossibleInDim(2)); + static FlatManifold<3,2> flat; + return flat; +} // explicit instantiations #include "tria.inst" diff --git a/tests/manifold/manifold_id_06.cc b/tests/manifold/manifold_id_06.cc index 9378c5f4fa..0a96f74a0a 100644 --- a/tests/manifold/manifold_id_06.cc +++ b/tests/manifold/manifold_id_06.cc @@ -58,7 +58,6 @@ void test(unsigned int ref=1) GridOut gridout; gridout.write_msh(tria, deallog.get_file_stream()); - } int main () diff --git a/tests/manifold/spherical_manifold_01.cc b/tests/manifold/spherical_manifold_01.cc index bbed8f87db..e645412452 100644 --- a/tests/manifold/spherical_manifold_01.cc +++ b/tests/manifold/spherical_manifold_01.cc @@ -32,7 +32,7 @@ void test(unsigned int ref=1) deallog << "Testing dim " << dim << ", spacedim " << spacedim << std::endl; - MINE::SphericalManifold manifold; + SphericalManifold manifold; Triangulation tria; GridGenerator::hyper_shell (tria, Point(), .3, .6, 12); diff --git a/tests/manifold/spherical_manifold_01.output b/tests/manifold/spherical_manifold_01.output index 59db6535ab..defe767cca 100644 --- a/tests/manifold/spherical_manifold_01.output +++ b/tests/manifold/spherical_manifold_01.output @@ -157,128 +157,128 @@ $NOD 26 0.00000 0.600000 0.00000 27 0.00000 0.00000 -0.600000 28 0.00000 0.00000 0.600000 -29 -0.264143 -0.109412 -0.0908716 +29 -0.266422 -0.0975173 -0.0975173 30 -0.259808 -0.259808 -0.259808 -31 0.264143 -0.109412 -0.0908716 -32 0.109412 -0.264143 -0.0908716 -33 0.127412 -0.0527760 -0.266422 +31 0.266422 -0.0975173 -0.0975173 +32 0.0975173 -0.266422 -0.0975173 +33 0.0975173 -0.0975173 -0.266422 34 0.259808 -0.259808 -0.259808 -35 -0.264143 0.109412 -0.0908716 -36 -0.109412 0.264143 -0.0908716 +35 -0.266422 0.0975173 -0.0975173 +36 -0.0975173 0.266422 -0.0975173 37 -0.259808 0.259808 -0.259808 -38 0.109412 0.264143 -0.0908716 +38 0.0975173 0.266422 -0.0975173 39 0.259808 0.259808 -0.259808 -40 -0.264143 -0.109412 0.0908716 -41 0.0527760 -0.127412 0.266422 +40 -0.266422 -0.0975173 0.0975173 +41 -0.0975173 -0.0975173 0.266422 42 -0.259808 -0.259808 0.259808 -43 0.127412 -0.0527760 0.266422 +43 0.0975173 -0.0975173 0.266422 44 0.259808 -0.259808 0.259808 45 -0.259808 0.259808 0.259808 -46 0.109412 0.264143 0.0908716 -47 0.127412 0.0527760 0.266422 +46 0.0975173 0.266422 0.0975173 +47 0.0975173 0.0975173 0.266422 48 0.259808 0.259808 0.259808 -49 -0.264143 0.109412 0.0908716 -50 -0.450000 5.51091e-17 2.75546e-17 -51 0.264143 0.109412 -0.0908716 -52 0.264143 -0.109412 0.0908716 -53 0.264143 0.109412 0.0908716 -54 0.450000 0.00000 2.75546e-17 -55 -0.109412 -0.264143 -0.0908716 -56 -0.109412 -0.264143 0.0908716 -57 0.109412 -0.264143 0.0908716 -58 -8.26637e-17 -0.450000 2.75546e-17 -59 -0.109412 0.264143 0.0908716 -60 2.75546e-17 0.450000 2.75546e-17 -61 0.0527760 -0.127412 -0.266422 -62 0.0527760 0.127412 -0.266422 -63 0.127412 0.0527760 -0.266422 -64 5.51091e-17 0.00000 -0.450000 -65 0.0527760 0.127412 0.266422 +49 -0.266422 0.0975173 0.0975173 +50 -0.450000 0.00000 0.00000 +51 0.266422 0.0975173 -0.0975173 +52 0.266422 -0.0975173 0.0975173 +53 0.266422 0.0975173 0.0975173 +54 0.450000 0.00000 0.00000 +55 -0.0975173 -0.266422 -0.0975173 +56 -0.0975173 -0.266422 0.0975173 +57 0.0975173 -0.266422 0.0975173 +58 0.00000 -0.450000 0.00000 +59 -0.0975173 0.266422 0.0975173 +60 0.00000 0.450000 0.00000 +61 -0.0975173 -0.0975173 -0.266422 +62 -0.0975173 0.0975173 -0.266422 +63 0.0975173 0.0975173 -0.266422 +64 0.00000 0.00000 -0.450000 +65 -0.0975173 0.0975173 0.266422 66 0.00000 0.00000 0.450000 -67 -0.528286 -0.218823 -0.181743 -68 0.528286 -0.218823 -0.181743 -69 0.218823 -0.528286 -0.181743 -70 0.254825 -0.105552 -0.532844 -71 -0.528286 0.218823 -0.181743 -72 -0.218823 0.528286 -0.181743 -73 0.218823 0.528286 -0.181743 -74 -0.528286 -0.218823 0.181743 -75 0.105552 -0.254825 0.532844 -76 0.254825 -0.105552 0.532844 -77 0.218823 0.528286 0.181743 -78 0.254825 0.105552 0.532844 -79 -0.528286 0.218823 0.181743 -80 0.528286 0.218823 -0.181743 -81 0.528286 -0.218823 0.181743 -82 0.528286 0.218823 0.181743 -83 -0.218823 -0.528286 -0.181743 -84 -0.218823 -0.528286 0.181743 -85 0.218823 -0.528286 0.181743 -86 -0.218823 0.528286 0.181743 -87 0.105552 -0.254825 -0.532844 -88 0.105552 0.254825 -0.532844 -89 0.254825 0.105552 -0.532844 -90 0.105552 0.254825 0.532844 -91 -0.396214 -0.164117 -0.136307 -92 0.229368 -5.61789e-17 -0.193366 -93 0.396214 -0.164117 -0.136307 -94 0.212132 -0.212132 1.83697e-17 -95 0.0877753 -0.211908 -0.193366 -96 0.191119 -0.0791640 -0.399633 -97 0.164117 -0.396214 -0.136307 -98 -0.396214 0.164117 -0.136307 -99 -0.212132 0.212132 1.83697e-17 -100 -0.164117 0.396214 -0.136307 -101 0.164117 0.396214 -0.136307 -102 0.211908 0.0877753 0.193366 -103 0.0791640 -0.191119 0.399633 -104 -0.396214 -0.164117 0.136307 -105 0.191119 -0.0791640 0.399633 -106 0.164117 0.396214 0.136307 -107 0.0877753 0.211908 0.193366 -108 0.191119 0.0791640 0.399633 -109 -0.396214 0.164117 0.136307 -110 0.396214 0.164117 -0.136307 -111 0.229368 -5.61789e-17 0.193366 -112 0.212132 0.212132 1.83697e-17 -113 0.396214 0.164117 0.136307 -114 0.396214 -0.164117 0.136307 -115 -0.212132 -0.212132 1.83697e-17 -116 0.0877753 -0.211908 0.193366 -117 -0.164117 -0.396214 0.136307 -118 0.164117 -0.396214 0.136307 -119 -0.164117 -0.396214 -0.136307 -120 -0.164117 0.396214 0.136307 -121 0.0791640 -0.191119 -0.399633 -122 0.211908 0.0877753 -0.193366 -123 0.0791640 0.191119 -0.399633 -124 0.0877753 0.211908 -0.193366 -125 0.191119 0.0791640 -0.399633 -126 0.0791640 0.191119 0.399633 -127 0.458736 -1.12358e-16 -0.386732 -128 0.424264 -0.424264 3.67394e-17 -129 0.175551 -0.423816 -0.386732 -130 -0.424264 0.424264 3.67394e-17 -131 0.423816 0.175551 0.386732 -132 0.175551 0.423816 0.386732 -133 0.458736 -1.12358e-16 0.386732 -134 0.424264 0.424264 3.67394e-17 -135 -0.424264 -0.424264 3.67394e-17 -136 0.175551 -0.423816 0.386732 -137 0.423816 0.175551 -0.386732 -138 0.175551 0.423816 -0.386732 -139 -0.318198 -0.318198 -1.72286e-16 -140 0.335064 0.0781260 0.290049 -141 0.131663 -0.317862 0.290049 -142 0.318198 -0.318198 -1.72286e-16 -143 0.344052 -8.42684e-17 0.290049 -144 0.131663 0.317862 0.290049 -145 0.318198 0.318198 2.75546e-17 -146 0.344052 -8.42684e-17 -0.290049 -147 0.131663 0.317862 -0.290049 -148 -0.318198 0.318198 -1.72286e-16 -149 0.335064 0.0781260 -0.290049 -150 0.131663 -0.317862 -0.290049 +67 -0.532844 -0.195035 -0.195035 +68 0.532844 -0.195035 -0.195035 +69 0.195035 -0.532844 -0.195035 +70 0.195035 -0.195035 -0.532844 +71 -0.532844 0.195035 -0.195035 +72 -0.195035 0.532844 -0.195035 +73 0.195035 0.532844 -0.195035 +74 -0.532844 -0.195035 0.195035 +75 -0.195035 -0.195035 0.532844 +76 0.195035 -0.195035 0.532844 +77 0.195035 0.532844 0.195035 +78 0.195035 0.195035 0.532844 +79 -0.532844 0.195035 0.195035 +80 0.532844 0.195035 -0.195035 +81 0.532844 -0.195035 0.195035 +82 0.532844 0.195035 0.195035 +83 -0.195035 -0.532844 -0.195035 +84 -0.195035 -0.532844 0.195035 +85 0.195035 -0.532844 0.195035 +86 -0.195035 0.532844 0.195035 +87 -0.195035 -0.195035 -0.532844 +88 -0.195035 0.195035 -0.532844 +89 0.195035 0.195035 -0.532844 +90 -0.195035 0.195035 0.532844 +91 -0.399633 -0.146276 -0.146276 +92 0.212132 0.00000 -0.212132 +93 0.399633 -0.146276 -0.146276 +94 0.212132 -0.212132 0.00000 +95 0.00000 -0.212132 -0.212132 +96 0.146276 -0.146276 -0.399633 +97 0.146276 -0.399633 -0.146276 +98 -0.399633 0.146276 -0.146276 +99 -0.212132 0.212132 0.00000 +100 -0.146276 0.399633 -0.146276 +101 0.146276 0.399633 -0.146276 +102 -0.212132 0.00000 0.212132 +103 -0.146276 -0.146276 0.399633 +104 -0.399633 -0.146276 0.146276 +105 0.146276 -0.146276 0.399633 +106 0.146276 0.399633 0.146276 +107 0.00000 0.212132 0.212132 +108 0.146276 0.146276 0.399633 +109 -0.399633 0.146276 0.146276 +110 0.399633 0.146276 -0.146276 +111 0.212132 0.00000 0.212132 +112 0.212132 0.212132 0.00000 +113 0.399633 0.146276 0.146276 +114 0.399633 -0.146276 0.146276 +115 -0.212132 -0.212132 0.00000 +116 0.00000 -0.212132 0.212132 +117 -0.146276 -0.399633 0.146276 +118 0.146276 -0.399633 0.146276 +119 -0.146276 -0.399633 -0.146276 +120 -0.146276 0.399633 0.146276 +121 -0.146276 -0.146276 -0.399633 +122 -0.212132 0.00000 -0.212132 +123 -0.146276 0.146276 -0.399633 +124 0.00000 0.212132 -0.212132 +125 0.146276 0.146276 -0.399633 +126 -0.146276 0.146276 0.399633 +127 0.424264 0.00000 -0.424264 +128 0.424264 -0.424264 0.00000 +129 0.00000 -0.424264 -0.424264 +130 -0.424264 0.424264 0.00000 +131 -0.424264 0.00000 0.424264 +132 0.00000 0.424264 0.424264 +133 0.424264 0.00000 0.424264 +134 0.424264 0.424264 0.00000 +135 -0.424264 -0.424264 0.00000 +136 0.00000 -0.424264 0.424264 +137 -0.424264 0.00000 -0.424264 +138 0.00000 0.424264 -0.424264 +139 -0.318198 -0.318198 0.00000 +140 -0.318198 0.00000 0.318198 +141 0.00000 -0.318198 0.318198 +142 0.318198 -0.318198 0.00000 +143 0.318198 0.00000 0.318198 +144 0.00000 0.318198 0.318198 +145 0.318198 0.318198 0.00000 +146 0.318198 0.00000 -0.318198 +147 0.00000 0.318198 -0.318198 +148 -0.318198 0.318198 0.00000 +149 -0.318198 0.00000 -0.318198 +150 0.00000 -0.318198 -0.318198 $ENDNOD $ELM 96 -- 2.39.5