From: Matthias Maier Date: Thu, 7 Aug 2014 11:15:32 +0000 (+0200) Subject: Take the chance to cleanup whitespace and indent X-Git-Tag: v8.2.0-rc1~220^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F64%2Fhead;p=dealii.git Take the chance to cleanup whitespace and indent --- diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h index b63c4aa981..b020a43e8d 100644 --- a/include/deal.II/grid/manifold.h +++ b/include/deal.II/grid/manifold.h @@ -36,7 +36,7 @@ template class Triangulation; * We collect here some helper functions used in the * Manifold classes. */ -namespace Manifolds +namespace Manifolds { /** * Given a hex iterator, construct a quadrature with the Laplace @@ -44,9 +44,9 @@ namespace Manifolds * centers and face centers, which can be called when creating * middle vertices in the manifold routines. */ - Quadrature<3> - get_default_quadrature(const TriaIterator >& hex); - + Quadrature<3> + get_default_quadrature(const TriaIterator > &hex); + /** * Given a general mesh iterator, construct a quadrature with the * Laplace weights or with uniform weights according the parameter @@ -55,8 +55,8 @@ namespace Manifolds * when creating new vertices in the manifold routines. */ template - Quadrature - get_default_quadrature(const OBJECT& obj, bool with_laplace = false); + Quadrature + get_default_quadrature(const OBJECT &obj, bool with_laplace = false); } @@ -83,7 +83,7 @@ namespace Manifolds * overload project_to_manifold() for simple situations. * * Should a finer control be necessary, then get_new_point() can be - * overloaded. + * overloaded. * * FlatManifold is the specialization from which StraigthBoundary is * derived, where the project_to_manifold() function is the identity. @@ -120,7 +120,7 @@ public: /** * 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 @@ -274,7 +274,7 @@ public: * argument. */ FlatManifold (const Point periodicity=Point(), - const double tolerance=1e-10); + const double tolerance=1e-10); /** * Let the new point be the average sum of surrounding vertices. @@ -300,8 +300,8 @@ public: */ 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 @@ -310,8 +310,8 @@ public: * similar (if not identical) to the one implemented in this class. */ virtual - Point project_to_manifold (const std::vector > &points, - const Point &candidate) const; + Point project_to_manifold (const std::vector > &points, + const Point &candidate) const; protected: /** * Tolerance. This tolerance is used to compute distances in double @@ -337,9 +337,9 @@ private: const Point periodicity; DeclException4(ExcPeriodicBox, int, Point, Point, double, - << "The component number " << arg1 << " of the point [ " << arg2 - << " ] is not in the interval [ " << -arg4 - << ", " << arg3[arg4] << "), bailing out."); + << "The component number " << arg1 << " of the point [ " << arg2 + << " ] is not in the interval [ " << -arg4 + << ", " << arg3[arg4] << "), bailing out."); }; @@ -379,7 +379,7 @@ private: * * Notice that the dimenisions #chartdim and #spacedim can be * arbitrary, as long as the transformation from $\mathcal{B}$ to - * $\mathcal{M}$ is invertible. + * $\mathcal{M}$ is invertible. * * @ingroup manifold * @@ -406,7 +406,7 @@ public: * 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. @@ -418,34 +418,34 @@ public: * 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 chartdim - * 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 chartdim dimensional 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: + virtual Point + get_new_point(const Quadrature &quad) const; + + /** + * Pull back the given point in spacedim to the Euclidean chartdim + * 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 chartdim dimensional 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. + * points in the chart coordinates system. */ const FlatManifold sub_manifold; }; @@ -498,61 +498,64 @@ get_new_point_on_hex (const Triangulation<3,3>::hex_iterator &) const; /*---Templated functions---*/ -namespace Manifolds { - +namespace Manifolds +{ + template - Quadrature - get_default_quadrature(const OBJECT& obj, bool with_laplace) + 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) + 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; + 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); } diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index fda040d710..8f54e8464f 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -24,7 +24,7 @@ DEAL_II_NAMESPACE_OPEN /** - * Manifold description for a spherical space coordinate system. + * Manifold description for a spherical space coordinate system. * * You can use this Manifold object to describe any sphere, circle, * hypersphere or hyperdisc in two or three dimensions, both as a @@ -47,13 +47,13 @@ DEAL_II_NAMESPACE_OPEN * be attached to cells containing the center of the coordinate * system. This point is a singular point of the coordinate * transformation, and there taking averages does not make any sense. - * + * * @ingroup manifold * * @author Luca Heltai, 2014 */ template -class SphericalManifold : public ManifoldChart +class SphericalManifold : public ManifoldChart { public: /** @@ -92,13 +92,13 @@ public: */ virtual Point get_new_point(const Quadrature &quad) const; - + /** * The center of the spherical coordinate system. */ const Point center; private: - + /** Helper function which returns the periodicity associated with this coordinate system, according to dim, chartdim, and spacedim. */ @@ -117,7 +117,7 @@ private: * This class was developed to be used in conjunction with the @p * cylinder or @p cylinder_shell functions of GridGenerator. This * function will throw an exception whenever spacedim is not equal to - * three. + * three. * * @ingroup manifold * @@ -134,8 +134,8 @@ public: * z-axis, respectively. The tolerance value is used to determine * if a point is on the axis. */ - CylindricalManifold (const unsigned int axis = 0, - const double tolerance = 1e-10); + CylindricalManifold (const unsigned int axis = 0, + const double tolerance = 1e-10); /** * Constructor. If constructed with this constructor, the manifold @@ -143,17 +143,17 @@ public: * #direction and goes through the given #point_on_axis. The * direction may be arbitrarily scaled, and the given point may be * any point on the axis. The tolerance value is used to determine - * if a point is on the axis. + * if a point is on the axis. */ CylindricalManifold (const Point &direction, - const Point &point_on_axis, - const double tolerance = 1e-10); + const Point &point_on_axis, + const double tolerance = 1e-10); - /** - * Compute new points on the CylindricalManifold. See the documentation - * of the base class for a detailed description of what this - * function does. - */ + /** + * Compute new points on the CylindricalManifold. See the documentation + * of the base class for a detailed description of what this + * function does. + */ virtual Point get_new_point(const Quadrature &quad) const; @@ -173,7 +173,7 @@ private: * Helper FlatManifold to compute temptative midpoints. */ FlatManifold flat_manifold; - + /** * Relative tolerance to measure zero distances. */ diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h index df114567ea..ed43e3f220 100644 --- a/include/deal.II/grid/tria.h +++ b/include/deal.II/grid/tria.h @@ -1244,7 +1244,7 @@ public: * boundary description has been explicitly * set using set_manifold(). */ - static const StraightBoundary straight_boundary; + static const StraightBoundary straight_boundary; /** * Declare some symbolic names @@ -1619,7 +1619,7 @@ public: virtual void create_notification (const Triangulation &tria); - }; + }; /** * A structure that is used as an @@ -1867,7 +1867,7 @@ public: void set_manifold (const types::manifold_id number, const Manifold &manifold_object); - + /** * Reset those parts of the triangulation with the given manifold_id * to use a FlatManifold object. This is the default state of a @@ -1897,7 +1897,7 @@ public: * @ingroup manifold * * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators" - */ + */ const Manifold &get_manifold (const types::manifold_id number) const; /** @@ -1912,7 +1912,7 @@ public: * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ std::vector get_boundary_indicators() const; - + /** * Returns a vector containing all manifold indicators assigned to * the objects of this Triangulation. Note, that each manifold @@ -1925,7 +1925,7 @@ public: * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators" */ std::vector get_manifold_ids() const; - + /** * Copy a triangulation. This operation is not cheap, so you should * be careful with using this. We do not implement this function as @@ -3375,7 +3375,7 @@ private: * a pointer. */ std::map *vertex_to_boundary_id_map_1d; - + /** * A map that relates the number of a boundary vertex to the @@ -3556,10 +3556,11 @@ Triangulation::save (Archive &ar, ar &check_for_distorted_cells; - if (dim == 1) { - ar &vertex_to_boundary_id_map_1d; - ar &vertex_to_manifold_id_map_1d; - } + if (dim == 1) + { + ar &vertex_to_boundary_id_map_1d; + ar &vertex_to_manifold_id_map_1d; + } } @@ -3592,10 +3593,11 @@ Triangulation::load (Archive &ar, "same setting with regard to reporting distorted " "cell as the one previously stored.")); - if (dim == 1) { - ar &vertex_to_boundary_id_map_1d; - ar &vertex_to_manifold_id_map_1d; - } + if (dim == 1) + { + ar &vertex_to_boundary_id_map_1d; + ar &vertex_to_manifold_id_map_1d; + } // trigger the create signal to indicate // that new content has been imported into @@ -3648,12 +3650,12 @@ 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; +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/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index 2f8c17fbdc..ace33fe1d1 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -1191,7 +1191,7 @@ public: * get_manifold() internally. */ const Boundary &get_boundary () const; - + /** * Return a constant reference to the manifold object used for this object. * @@ -1269,8 +1269,8 @@ public: /** * @} - */ - + */ + /** * @name User data @@ -2106,7 +2106,7 @@ public: /** * Return the manifold indicator of this - * object. + * object. * * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators" */ @@ -2285,7 +2285,7 @@ public: * * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators" */ - + void set_all_manifold_ids (const types::manifold_id); /** diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index d6cda3e65f..842bd16337 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -1912,21 +1912,21 @@ TriaAccessor::at_boundary () const template const Boundary & TriaAccessor::get_boundary () const -{ +{ Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - - // Get the default (manifold_id) + + // Get the default (manifold_id) const types::manifold_id mi = this->objects().manifold_id[this->present_index]; - - // In case this is not valid, check - // the boundary id, after having - // casted it to a manifold id - if(mi == numbers::invalid_manifold_id) + + // In case this is not valid, check + // the boundary id, after having + // casted it to a manifold id + if (mi == numbers::invalid_manifold_id) return this->tria->get_boundary(structdim < dim ? - this->objects().boundary_or_material_id[this->present_index].boundary_id: - dim < spacedim ? - this->objects().boundary_or_material_id[this->present_index].material_id: - numbers::invalid_manifold_id); + this->objects().boundary_or_material_id[this->present_index].boundary_id: + dim < spacedim ? + this->objects().boundary_or_material_id[this->present_index].material_id: + numbers::invalid_manifold_id); else return this->tria->get_boundary(mi); } @@ -1935,21 +1935,21 @@ TriaAccessor::get_boundary () const template const Manifold & TriaAccessor::get_manifold () const -{ +{ Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - - // Get the default (manifold_id) + + // Get the default (manifold_id) const types::manifold_id mi = this->objects().manifold_id[this->present_index]; - - // In case this is not valid, check - // the boundary id, after having - // casted it to a manifold id - if(mi == numbers::invalid_manifold_id) + + // In case this is not valid, check + // the boundary id, after having + // casted it to a manifold id + if (mi == numbers::invalid_manifold_id) return this->tria->get_manifold(structdim < dim ? - this->objects().boundary_or_material_id[this->present_index].boundary_id: - dim < spacedim ? - this->objects().boundary_or_material_id[this->present_index].material_id: - numbers::invalid_manifold_id); + this->objects().boundary_or_material_id[this->present_index].boundary_id: + dim < spacedim ? + this->objects().boundary_or_material_id[this->present_index].material_id: + numbers::invalid_manifold_id); else return this->tria->get_manifold(mi); } @@ -1991,13 +1991,13 @@ set_all_manifold_ids (const types::manifold_id manifold_ind) const switch (structdim) { case 1: - if(dim == 1) - { - (*this->tria->vertex_to_manifold_id_map_1d) - [vertex_index(0)] = manifold_ind; - (*this->tria->vertex_to_manifold_id_map_1d) - [vertex_index(1)] = manifold_ind; - } + if (dim == 1) + { + (*this->tria->vertex_to_manifold_id_map_1d) + [vertex_index(0)] = manifold_ind; + (*this->tria->vertex_to_manifold_id_map_1d) + [vertex_index(1)] = manifold_ind; + } break; case 2: @@ -2376,8 +2376,8 @@ inline types::manifold_id TriaAccessor<0, 1, spacedim>::manifold_id () const { - if( tria->vertex_to_manifold_id_map_1d->find (this->vertex_index()) - != tria->vertex_to_manifold_id_map_1d->end()) + if ( tria->vertex_to_manifold_id_map_1d->find (this->vertex_index()) + != tria->vertex_to_manifold_id_map_1d->end()) return (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()]; else return numbers::invalid_manifold_id; diff --git a/include/deal.II/grid/tria_objects.h b/include/deal.II/grid/tria_objects.h index db5c01d685..e7b0350875 100644 --- a/include/deal.II/grid/tria_objects.h +++ b/include/deal.II/grid/tria_objects.h @@ -206,7 +206,7 @@ namespace internal * Store manifold ids. This field * stores the manifold id of each object, which * is a number between 0 and - * numbers::invalid_manifold_id-1. + * numbers::invalid_manifold_id-1. */ std::vector manifold_id; diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index b16ab5765f..0b728dacb3 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -24,47 +24,48 @@ DEAL_II_NAMESPACE_OPEN -namespace Manifolds { +namespace Manifolds +{ - Quadrature<3> - get_default_quadrature(const TriaIterator >& obj) + Quadrature<3> + get_default_quadrature(const TriaIterator > &obj) { - std::vector > sp; - std::vector wp; - - const int dim = 3; - - 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); - 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)->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)->get_manifold().get_new_point_on_face(obj->face(i))); - wp[j] = 1.0/12.0; - } - return Quadrature<3>(sp,wp); + std::vector > sp; + std::vector wp; + + const int dim = 3; + + 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); + 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)->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)->get_manifold().get_new_point_on_face(obj->face(i))); + wp[j] = 1.0/12.0; + } + return Quadrature<3>(sp,wp); } } @@ -84,7 +85,7 @@ template Point Manifold:: project_to_manifold (const std::vector > &, - const Point &) const + const Point &) const { Assert (false, ExcPureFunctionCalled()); return Point(); @@ -103,12 +104,12 @@ get_new_point (const Quadrature &quad) const #ifdef DEBUG double sum=0; - for(unsigned int i=0; i:: get_new_point_on_line (const typename Triangulation::line_iterator &line) const { return get_new_point - (get_default_quadrature::line_iterator, - spacedim>(line, false)); + (get_default_quadrature::line_iterator, + spacedim>(line, false)); } @@ -133,8 +134,8 @@ Manifold:: get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const { return get_new_point - (get_default_quadrature::quad_iterator, - spacedim>(quad, false)); + (get_default_quadrature::quad_iterator, + spacedim>(quad, false)); } @@ -248,7 +249,8 @@ get_new_point_on_hex (const typename Triangulation::hex_iterator template <> Point<3> Manifold<3,3>:: -get_new_point_on_hex (const typename Triangulation<3, 3>::hex_iterator &hex) const{ +get_new_point_on_hex (const typename Triangulation<3, 3>::hex_iterator &hex) const +{ return get_new_point(get_default_quadrature(hex)); } @@ -257,8 +259,8 @@ get_new_point_on_hex (const typename Triangulation<3, 3>::hex_iterator &hex) con template -FlatManifold::FlatManifold (const Point periodicity, - const double tolerance) : +FlatManifold::FlatManifold (const Point periodicity, + const double tolerance) : tolerance(tolerance), periodicity(periodicity) {} @@ -273,48 +275,51 @@ get_new_point (const Quadrature &quad) const #ifdef DEBUG double sum=0; - for(unsigned int i=0; i p; Point dp; Point minP = periodicity; const bool check_period = (periodicity.norm() > tolerance); - if(check_period) - for(unsigned int i=0; i 0) - Assert( (surrounding_points[i][d] < periodicity[d]+tolerance) || - (surrounding_points[i][d] >= -tolerance), - ExcPeriodicBox(d, surrounding_points[i], periodicity, tolerance)); - } - - 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 ); + if (check_period) + for (unsigned int i=0; i 0) + Assert( (surrounding_points[i][d] < periodicity[d]+tolerance) || + (surrounding_points[i][d] >= -tolerance), + ExcPeriodicBox(d, surrounding_points[i], periodicity, tolerance)); + } + + 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]; } - 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]); + 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 +Point +FlatManifold::project_to_manifold (const std::vector > &vertices, + const Point &candidate) const { return candidate; } @@ -340,13 +345,13 @@ 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 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 bd479a1255..d2781df4d4 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -35,7 +35,8 @@ SphericalManifold::SphericalManifold(const Point center) template Point -SphericalManifold::get_periodicity() { +SphericalManifold::get_periodicity() +{ Point periodicity; periodicity[spacedim-1] = 2*numbers::PI; // theta and phi period. return periodicity; @@ -44,85 +45,92 @@ SphericalManifold::get_periodicity() { template Point -SphericalManifold::get_new_point(const Quadrature &quad) const { - if(spacedim == 2) +SphericalManifold::get_new_point(const Quadrature &quad) const +{ + if (spacedim == 2) return ManifoldChart::get_new_point(quad); - else { - double rho_average = 0; - Point mid_point; - for(unsigned int i=0; i mid_point; + for (unsigned int i=0; i R = mid_point-center; + // Scale it to have radius rho_average + R *= rho_average/R.norm(); + // And return it. + return center+R; } - // Project the mid_pont back to the right location - Point R = mid_point-center; - // Scale it to have radius rho_average - R *= rho_average/R.norm(); - // And return it. - return center+R; - } } template Point -SphericalManifold::push_forward(const Point &spherical_point) const { - Assert(spherical_point[0] >=0.0, - ExcMessage("Negative radius for given 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]; - + Point p; - if(rho > 1e-10) - switch(spacedim) { - case 2: - p[0] = rho*cos(theta); - p[1] = rho*sin(theta); - break; - case 3: + if (rho > 1e-10) + switch (spacedim) + { + case 2: + p[0] = rho*cos(theta); + p[1] = rho*sin(theta); + break; + case 3: { - const double &phi= spherical_point[2]; - p[0] = rho*sin(theta)*cos(phi); - p[1] = rho*sin(theta)*sin(phi); - p[2] = rho*cos(theta); + const double &phi= spherical_point[2]; + p[0] = rho*sin(theta)*cos(phi); + p[1] = rho*sin(theta)*sin(phi); + p[2] = rho*cos(theta); } break; - default: - Assert(false, ExcInternalError()); - } + default: + Assert(false, ExcInternalError()); + } return p+center; } template Point -SphericalManifold::pull_back(const Point &space_point) const { +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]; - + Point p; p[0] = rho; - - switch(spacedim) { - case 2: - p[1] = atan2(y,x); - if(p[1] < 0) - p[1] += 2*numbers::PI; - break; - case 3: + + switch (spacedim) + { + case 2: + p[1] = atan2(y,x); + if (p[1] < 0) + p[1] += 2*numbers::PI; + break; + case 3: { const double z = R[2]; p[2] = atan2(y,x); // phi - if(p[2] < 0) - p[2] += 2*numbers::PI; // phi is periodic + if (p[2] < 0) + p[2] += 2*numbers::PI; // phi is periodic p[1] = atan2(sqrt(x*x+y*y),z); // theta } break; - default: - Assert(false, ExcInternalError()); - } + default: + Assert(false, ExcInternalError()); + } return p; } @@ -132,8 +140,8 @@ SphericalManifold::pull_back(const Point &space_point) c // ============================================================ template -CylindricalManifold::CylindricalManifold(const unsigned int axis, - const double tolerance) : +CylindricalManifold::CylindricalManifold(const unsigned int axis, + const double tolerance) : direction (Point::unit_vector(axis)), point_on_axis (Point()), tolerance(tolerance) @@ -144,8 +152,8 @@ CylindricalManifold::CylindricalManifold(const unsigned int axis, template CylindricalManifold::CylindricalManifold(const Point &direction, - const Point &point_on_axis, - const double tolerance) : + const Point &point_on_axis, + const double tolerance) : direction (direction), point_on_axis (point_on_axis), tolerance(tolerance) @@ -164,24 +172,24 @@ get_new_point (const Quadrature &quad) const const std::vector > &surrounding_points = quad.get_points(); const std::vector &weights = quad.get_weights(); - // compute a proposed new point + // compute a proposed new point Point middle = flat_manifold.get_new_point(quad); double radius = 0; Point on_plane; - - for(unsigned int i=0; i vector_from_axis = (middle-point_on_axis) - - ((middle-point_on_axis) * direction) * direction; + ((middle-point_on_axis) * direction) * direction; // scale it to the desired length and put everything back together, // unless we have a point on the axis @@ -190,10 +198,10 @@ get_new_point (const Quadrature &quad) const else return (vector_from_axis / vector_from_axis.norm() * radius + - ((middle-point_on_axis) * direction) * direction + - point_on_axis); + ((middle-point_on_axis) * direction) * direction + + point_on_axis); } - + // explicit instantiations #include "manifold_lib.inst" diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 8169c5952a..03380582b9 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -1488,7 +1488,7 @@ namespace internal cells[cell].vertices[1])); next_free_line->set_used_flag (); next_free_line->set_material_id (cells[cell].material_id); - next_free_line->set_manifold_id (cells[cell].manifold_id); + next_free_line->set_manifold_id (cells[cell].manifold_id); next_free_line->clear_user_data (); next_free_line->set_subdomain_id (0); @@ -1576,24 +1576,24 @@ namespace internal // finally set the // vertex_to_boundary_id_map_1d - // and vertex_to_manifold_id_map_1d + // and vertex_to_manifold_id_map_1d // maps triangulation.vertex_to_boundary_id_map_1d->clear(); - triangulation.vertex_to_manifold_id_map_1d->clear(); + triangulation.vertex_to_manifold_id_map_1d->clear(); for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end(); ++cell) for (unsigned int f=0; f::faces_per_cell; ++f) - { - (*triangulation - .vertex_to_manifold_id_map_1d)[cell->face(f)->vertex_index()] - = numbers::invalid_manifold_id; - - if (cell->at_boundary(f)) - (*triangulation - .vertex_to_boundary_id_map_1d)[cell->face(f)->vertex_index()] - = f; - } + { + (*triangulation + .vertex_to_manifold_id_map_1d)[cell->face(f)->vertex_index()] + = numbers::invalid_manifold_id; + + if (cell->at_boundary(f)) + (*triangulation + .vertex_to_boundary_id_map_1d)[cell->face(f)->vertex_index()] + = f; + } } @@ -1770,7 +1770,7 @@ namespace internal cell->set_used_flag (); cell->set_material_id (cells[c].material_id); - cell->set_manifold_id (cells[c].manifold_id); + cell->set_manifold_id (cells[c].manifold_id); cell->clear_user_data (); cell->set_subdomain_id (0); @@ -1803,7 +1803,7 @@ namespace internal else // interior line -> numbers::internal_face_boundary_id line->set_boundary_indicator (numbers::internal_face_boundary_id); - line->set_manifold_id(numbers::flat_manifold_id); + line->set_manifold_id(numbers::flat_manifold_id); } // set boundary indicators where @@ -1847,7 +1847,7 @@ namespace internal ExcInteriorLineCantBeBoundary()); line->set_boundary_indicator (boundary_line->boundary_id); - line->set_manifold_id (boundary_line->manifold_id); + line->set_manifold_id (boundary_line->manifold_id); } @@ -2424,7 +2424,7 @@ namespace internal cell->set_used_flag (); cell->set_material_id (cells[c].material_id); - cell->set_manifold_id (cells[c].manifold_id); + cell->set_manifold_id (cells[c].manifold_id); cell->clear_user_flag (); cell->clear_user_data (); cell->set_subdomain_id (0); @@ -2526,10 +2526,10 @@ namespace internal else // interior quad -> numbers::internal_face_boundary_id quad->set_boundary_indicator (numbers::internal_face_boundary_id); - // Manifold ids are set - // independently of where - // they are - quad->set_manifold_id(numbers::flat_manifold_id); + // Manifold ids are set + // independently of where + // they are + quad->set_manifold_id(numbers::flat_manifold_id); } ///////////////////////////////////////// @@ -2540,12 +2540,12 @@ namespace internal // for this: first mark all lines as interior. use this loop // to also set all manifold ids of all lines for (typename Triangulation::line_iterator - line=triangulation.begin_line(); line!=triangulation.end_line(); ++line) - { - line->set_boundary_indicator (numbers::internal_face_boundary_id); - line->set_manifold_id(numbers::flat_manifold_id); - } - + line=triangulation.begin_line(); line!=triangulation.end_line(); ++line) + { + line->set_boundary_indicator (numbers::internal_face_boundary_id); + line->set_manifold_id(numbers::flat_manifold_id); + } + // next reset all lines bounding // boundary quads as on the // boundary also. note that since @@ -2618,8 +2618,8 @@ namespace internal "if they carry the same boundary indicator.")); line->set_boundary_indicator (boundary_line->boundary_id); - // Set manifold id if given - line->set_manifold_id(boundary_line->manifold_id); + // Set manifold id if given + line->set_manifold_id(boundary_line->manifold_id); } @@ -2752,7 +2752,7 @@ namespace internal "if they carry the same boundary indicator.")); quad->set_boundary_indicator (boundary_quad->boundary_id); - quad->set_manifold_id (boundary_quad->manifold_id); + quad->set_manifold_id (boundary_quad->manifold_id); } @@ -3607,7 +3607,7 @@ namespace internal switch_1->set_line_orientation(2, switch_2->line_orientation(2)); switch_1->set_line_orientation(3, switch_2->line_orientation(3)); switch_1->set_boundary_indicator(switch_2->boundary_indicator()); - switch_1->set_manifold_id(switch_2->manifold_id()); + switch_1->set_manifold_id(switch_2->manifold_id()); switch_1->set_user_index(switch_2->user_index()); if (switch_2->user_flag_set()) switch_1->set_user_flag(); @@ -3623,7 +3623,7 @@ namespace internal switch_2->set_line_orientation(2, switch_1_line_orientations[2]); switch_2->set_line_orientation(3, switch_1_line_orientations[3]); switch_2->set_boundary_indicator(switch_1_boundary_indicator); - switch_2->set_manifold_id(switch_1->manifold_id()); + switch_2->set_manifold_id(switch_1->manifold_id()); switch_2->set_user_index(switch_1_user_index); if (switch_1_user_flag) switch_2->set_user_flag(); @@ -3955,11 +3955,11 @@ namespace internal if (dim == spacedim) { // triangulation.vertices[next_unused_vertex] = new_point; - triangulation.vertices[next_unused_vertex] = - cell->get_manifold().get_new_point - (Manifolds::get_default_quadrature - ::active_cell_iterator, - spacedim>(cell)); + triangulation.vertices[next_unused_vertex] = + 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 @@ -4004,19 +4004,19 @@ namespace internal // points on this face and on the opposite face, // as returned by the underlying manifold // object. - { - 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); - } - } + { + 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 { @@ -4037,8 +4037,8 @@ namespace internal // new vertex is placed on the surface according to // the information stored in the boundary class - const Manifold &manifold = cell->get_manifold(); - + const Manifold &manifold = cell->get_manifold(); + triangulation.vertices[next_unused_vertex] = manifold.get_new_point_on_quad (cell); } @@ -4133,7 +4133,7 @@ namespace internal new_lines[l]->clear_children(); // interior line new_lines[l]->set_boundary_indicator(numbers::internal_face_boundary_id); - new_lines[l]->set_manifold_id(cell->manifold_id()); + new_lines[l]->set_manifold_id(cell->manifold_id()); } // Now add the four (two) @@ -4255,7 +4255,7 @@ namespace internal // inherit material // properties subcells[i]->set_material_id (cell->material_id()); - subcells[i]->set_manifold_id (cell->manifold_id()); + subcells[i]->set_manifold_id (cell->manifold_id()); subcells[i]->set_subdomain_id (subdomainid); if (i%2==0) @@ -4403,19 +4403,19 @@ namespace internal Assert (next_unused_vertex < triangulation.vertices.size(), ExcTooFewVerticesAllocated()); - // Now we always ask the manifold where to put the - // new point. The get_manifold function will return - // a flat boundary if invalid_manifold_id is set, - // behaving as before in the flat case. Backward - // compatibility requires us to use the material id - // of the cell in the codimension one case, however - // we only do this if the manifold_id is the invalid - // one, otherwise use the material id. This is done - // internally in the get_manifold() function. - triangulation.vertices[next_unused_vertex] = - cell->get_manifold().get_new_point_on_line(cell); - triangulation.vertices_used[next_unused_vertex] = true; - + // Now we always ask the manifold where to put the + // new point. The get_manifold function will return + // a flat boundary if invalid_manifold_id is set, + // behaving as before in the flat case. Backward + // compatibility requires us to use the material id + // of the cell in the codimension one case, however + // we only do this if the manifold_id is the invalid + // one, otherwise use the material id. This is done + // internally in the get_manifold() function. + triangulation.vertices[next_unused_vertex] = + cell->get_manifold().get_new_point_on_line(cell); + triangulation.vertices_used[next_unused_vertex] = true; + // search for next two unused cell (++ takes care of // the end of the vector) typename Triangulation::raw_cell_iterator @@ -4442,16 +4442,16 @@ namespace internal ::TriaObject<1> (cell->vertex_index(0), next_unused_vertex)); first_child->set_material_id (cell->material_id()); - first_child->set_manifold_id (cell->manifold_id()); + first_child->set_manifold_id (cell->manifold_id()); first_child->set_subdomain_id (subdomainid); first_child->set_direction_flag (cell->direction_flag()); - + first_child->set_parent (cell->index ()); - // Set manifold id of the right face. Only do this - // on the first child. - first_child->face(1)->set_manifold_id(cell->manifold_id()); - + // Set manifold id of the right face. Only do this + // on the first child. + first_child->face(1)->set_manifold_id(cell->manifold_id()); + // reset neighborship info (refer to // internal::Triangulation::TriaLevel<0> for // details) @@ -4493,7 +4493,7 @@ namespace internal cell->vertex_index(1))); second_child->set_neighbor (0, first_child); second_child->set_material_id (cell->material_id()); - second_child->set_manifold_id (cell->manifold_id()); + second_child->set_manifold_id (cell->manifold_id()); second_child->set_subdomain_id (subdomainid); second_child->set_direction_flag (cell->direction_flag()); @@ -4746,8 +4746,8 @@ namespace internal // boundary lines differently; for interior // lines we can compute the midpoint as the mean // of the two vertices: if (line->at_boundary()) - triangulation.vertices[next_unused_vertex] - = line->get_manifold().get_new_point_on_line (line); + triangulation.vertices[next_unused_vertex] + = line->get_manifold().get_new_point_on_line (line); } else // however, if spacedim>dim, we always have to ask @@ -4755,13 +4755,13 @@ namespace internal // same object of the cell (which was stored in // line->user_index() before) unless a manifold_id // has been set on this very line. - if(line->manifold_id() == numbers::invalid_manifold_id) - triangulation.vertices[next_unused_vertex] - = triangulation.get_manifold(line->user_index()).get_new_point_on_line (line); - else - triangulation.vertices[next_unused_vertex] - = line->get_manifold().get_new_point_on_line (line); - + if (line->manifold_id() == numbers::invalid_manifold_id) + triangulation.vertices[next_unused_vertex] + = triangulation.get_manifold(line->user_index()).get_new_point_on_line (line); + else + triangulation.vertices[next_unused_vertex] + = line->get_manifold().get_new_point_on_line (line); + // now that we created the right point, make up the // two child lines. To this end, find a pair of // unused lines @@ -5204,8 +5204,8 @@ namespace internal ExcTooFewVerticesAllocated()); triangulation.vertices_used[next_unused_vertex] = true; - triangulation.vertices[next_unused_vertex] - = line->get_manifold().get_new_point_on_line (line); + triangulation.vertices[next_unused_vertex] + = line->get_manifold().get_new_point_on_line (line); // now that we created the right point, make up the // two child lines (++ takes care of the end of the @@ -5248,8 +5248,8 @@ namespace internal children[0]->set_boundary_indicator (line->boundary_indicator()); children[1]->set_boundary_indicator (line->boundary_indicator()); - - children[0]->set_manifold_id (line->manifold_id()); + + children[0]->set_manifold_id (line->manifold_id()); children[1]->set_manifold_id (line->manifold_id()); // finally clear flag @@ -5362,7 +5362,7 @@ namespace internal new_line->clear_user_data(); new_line->clear_children(); new_line->set_boundary_indicator(quad->boundary_indicator()); - new_line->set_manifold_id(quad->manifold_id()); + new_line->set_manifold_id(quad->manifold_id()); // child 0 and 1 of a line are switched if the // line orientation is false. set up a miniature @@ -5423,7 +5423,7 @@ namespace internal new_quads[i]->clear_user_data(); new_quads[i]->clear_children(); new_quads[i]->set_boundary_indicator (quad->boundary_indicator()); - new_quads[i]->set_manifold_id (quad->manifold_id()); + new_quads[i]->set_manifold_id (quad->manifold_id()); // set all line orientations to true, change // this after the loop, as we have to consider // different lines for each child @@ -5524,7 +5524,7 @@ namespace internal new_child[i]->set(internal::Triangulation::TriaObject<1>(old_child[i]->vertex_index(0), old_child[i]->vertex_index(1))); new_child[i]->set_boundary_indicator(old_child[i]->boundary_indicator()); - new_child[i]->set_manifold_id(old_child[i]->manifold_id()); + new_child[i]->set_manifold_id(old_child[i]->manifold_id()); new_child[i]->set_user_index(old_child[i]->user_index()); if (old_child[i]->user_flag_set()) new_child[i]->set_user_flag(); @@ -5627,7 +5627,7 @@ namespace internal switch_1->set_line_orientation(2, switch_2->line_orientation(2)); switch_1->set_line_orientation(3, switch_2->line_orientation(3)); switch_1->set_boundary_indicator(switch_2->boundary_indicator()); - switch_1->set_manifold_id(switch_2->manifold_id()); + switch_1->set_manifold_id(switch_2->manifold_id()); switch_1->set_user_index(switch_2->user_index()); if (switch_2->user_flag_set()) switch_1->set_user_flag(); @@ -5650,7 +5650,7 @@ namespace internal switch_2->set_line_orientation(2, switch_1_line_orientations[2]); switch_2->set_line_orientation(3, switch_1_line_orientations[3]); switch_2->set_boundary_indicator(switch_1_boundary_indicator); - switch_2->set_manifold_id(switch_1->manifold_id()); + switch_2->set_manifold_id(switch_1->manifold_id()); switch_2->set_user_index(switch_1_user_index); if (switch_1_user_flag) switch_2->set_user_flag(); @@ -5798,40 +5798,40 @@ namespace internal // set the middle vertex // appropriately - if (quad->at_boundary() || - (quad->manifold_id() != numbers::invalid_manifold_id) ) + if (quad->at_boundary() || + (quad->manifold_id() != numbers::invalid_manifold_id) ) triangulation.vertices[next_unused_vertex] - = quad->get_manifold().get_new_point_on_quad (quad); + = 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->get_manifold().get_new_point - (Manifolds::get_default_quadrature - ::quad_iterator, - spacedim>(quad, true)); - } - 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->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 // the four lines interior to the quad (++ takes // care of the end of the vector) @@ -5865,7 +5865,7 @@ namespace internal // *--2--* // // the lines are numbered as follows: - // + // // *--*--* // | 1 | // *2-*-3* @@ -5896,13 +5896,13 @@ namespace internal new_lines[i]->clear_user_data(); new_lines[i]->clear_children(); new_lines[i]->set_boundary_indicator(quad->boundary_indicator()); - new_lines[i]->set_manifold_id(quad->manifold_id()); + new_lines[i]->set_manifold_id(quad->manifold_id()); } // now for the quads. again, first collect some // data about the indices of the lines, with the // following numbering: - // + // // .-6-.-7-. // 1 9 3 // .-10.11-. @@ -5996,7 +5996,7 @@ namespace internal new_quads[i]->clear_user_data(); new_quads[i]->clear_children(); new_quads[i]->set_boundary_indicator (quad->boundary_indicator()); - new_quads[i]->set_manifold_id (quad->manifold_id()); + new_quads[i]->set_manifold_id (quad->manifold_id()); // set all line orientations to true, change // this after the loop, as we have to consider // different lines for each child @@ -6129,7 +6129,7 @@ namespace internal // interior quad new_quads[i]->set_boundary_indicator (numbers::internal_face_boundary_id); // they inherit geometry description of the hex they belong to - new_quads[i]->set_manifold_id (hex->manifold_id()); + new_quads[i]->set_manifold_id (hex->manifold_id()); // set all line orientation flags to true by // default, change this afterwards, if necessary for (unsigned int j=0; j::lines_per_face; ++j) @@ -6160,7 +6160,7 @@ namespace internal // inherit material // properties new_hexes[i]->set_material_id (hex->material_id()); - new_hexes[i]->set_manifold_id (hex->manifold_id()); + new_hexes[i]->set_manifold_id (hex->manifold_id()); new_hexes[i]->set_subdomain_id (subdomainid); if (i%2) @@ -6393,7 +6393,7 @@ namespace internal // the quads are numbered as follows: // // planes in the interior of the old hex: - // + // // * // /| // / | x @@ -6407,7 +6407,7 @@ namespace internal // * // // children of the faces of the old hex - // + // // *---*---* *---*---* // /| | | / / /| // / | | | / 9 / 10/ | @@ -6593,7 +6593,7 @@ namespace internal // the quads are numbered as follows: // // planes in the interior of the old hex: - // + // // * // /| // / | x @@ -6607,7 +6607,7 @@ namespace internal // * // // children of the faces of the old hex - // + // // *-------* *-------* // /| | / 10 /| // * | | *-------* | @@ -6795,7 +6795,7 @@ namespace internal // the quads are numbered as follows: // // planes in the interior of the old hex: - // + // // * // /| // / | x @@ -6809,7 +6809,7 @@ namespace internal // * // // children of the faces of the old hex - // + // // *---*---* *-------* // /| 8 | / /| // / | | / 10 / | @@ -7115,7 +7115,7 @@ namespace internal // the quads are numbered as follows: // // planes in the interior of the old hex: - // + // // * // /| // * | x @@ -7129,7 +7129,7 @@ namespace internal // * // // children of the faces of the old hex - // + // // *---*---* *---*---* // /| | | /18 / 19/| // * |10 | 11| /---/---* | @@ -7457,7 +7457,7 @@ namespace internal // the quads are numbered as follows: // // planes in the interior of the old hex: - // + // // * // /| // / | x @@ -7805,7 +7805,7 @@ namespace internal // the quads are numbered as follows: // // planes in the interior of the old hex: - // + // // * // /| // / | x @@ -7936,34 +7936,15 @@ namespace internal // function, and like the new 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] = - hex->get_manifold().get_new_point_on_hex(hex); + triangulation.vertices[next_unused_vertex] = + hex->get_manifold().get_new_point_on_hex(hex); // set the data of the six lines. first collect // the indices of the seven vertices (consider // the two planes to be crossed to form the // planes cutting the hex in two vertically and // horizontally) - // + // // *--3--* *--5--* // / / / | | | // 0--6--1 0--6--1 @@ -8309,7 +8290,7 @@ namespace internal // the quads are numbered as follows: // // planes in the interior of the old hex: - // + // // * // /| // * | @@ -8904,9 +8885,9 @@ Triangulation (const MeshSmoothing smooth_grid, if (dim == 1) { vertex_to_boundary_id_map_1d - = new std::map(); + = new std::map(); vertex_to_manifold_id_map_1d - = new std::map(); + = new std::map(); } // connect the any_change signal to the other signals @@ -8994,7 +8975,7 @@ Triangulation::set_manifold (const types::manifold_id m_number, const Manifold &manifold_object) { Assert(m_number < numbers::invalid_manifold_id, - ExcIndexRange(m_number,0,numbers::invalid_manifold_id)); + ExcIndexRange(m_number,0,numbers::invalid_manifold_id)); manifold[m_number] = &manifold_object; } @@ -9013,7 +8994,7 @@ void Triangulation::set_manifold (const types::manifold_id m_number) { Assert(m_number < numbers::invalid_manifold_id, - ExcIndexRange(m_number,0,numbers::invalid_manifold_id)); + ExcIndexRange(m_number,0,numbers::invalid_manifold_id)); //delete the entry located at number. manifold.erase(m_number); @@ -9025,11 +9006,11 @@ template const Boundary & Triangulation::get_boundary (const types::manifold_id m_number) const { - const Boundary * man = + const Boundary *man = dynamic_cast *>(&get_manifold(m_number)); Assert(man != NULL, - ExcMessage("You tried to get a Boundary, but I only have a Manifold.")); - + ExcMessage("You tried to get a Boundary, but I only have a Manifold.")); + return *man; } @@ -9038,8 +9019,8 @@ template const Manifold & Triangulation::get_manifold (const types::manifold_id m_number) const { - //look, if there is a manifold stored at - //manifold_id number. + //look, if there is a manifold stored at + //manifold_id number. typename std::map, Triangulation > >::const_iterator it = manifold.find(m_number); @@ -9095,13 +9076,14 @@ Triangulation::get_manifold_ids () const { std::set m_ids; active_cell_iterator cell=begin_active(); - for (; cell!=end(); ++cell) { - m_ids.insert(cell->manifold_id()); - if(dim>1) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->at_boundary(face)) - m_ids.insert(cell->face(face)->manifold_id()); - } + for (; cell!=end(); ++cell) + { + m_ids.insert(cell->manifold_id()); + if (dim>1) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->at_boundary(face)) + m_ids.insert(cell->face(face)->manifold_id()); + } std::vector manifold_indicators(m_ids.begin(), m_ids.end()); return manifold_indicators; } @@ -9132,8 +9114,8 @@ copy_triangulation (const Triangulation &old_tria) faces = new internal::Triangulation::TriaFaces(*old_tria.faces); typename std::map , Triangulation > >::const_iterator - bdry_iterator = old_tria.manifold.begin(); + SmartPointer , Triangulation > >::const_iterator + bdry_iterator = old_tria.manifold.begin(); for (; bdry_iterator != old_tria.manifold.end() ; bdry_iterator++) manifold[bdry_iterator->first] = bdry_iterator->second; @@ -11075,7 +11057,7 @@ unsigned int Triangulation::n_active_cells (const unsigned int le template bool Triangulation::has_hanging_nodes () const { - for (unsigned int lvl = 0; lvl::prepare_coarsening_and_refinement () // coarsen flags on level 1 to avoid level 0 cells being created // by coarsening. As coarsen flags will never be added, this can // be done once and for all before the actual loop starts. - // + // // 1/ do not coarsen a cell if 'most of the neighbors' will be // refined after the step. This is to prevent occurrence of // unrefined islands. - // + // // 2/ eliminate refined islands in the interior and at the // boundary. since they don't do much harm besides increasing the // number of degrees of freedom, doing this has a rather low // priority. - // + // // 3/ limit the level difference of neighboring cells at each // vertex. - // + // // 4/ eliminate unrefined islands. this has higher priority since // this diminishes the approximation properties not only of the // unrefined island, but also of the surrounding patch. - // + // // 5/ ensure patch level 1. Then the triangulation consists of // patches, i.e. of cells that are refined once. It follows that // if at least one of the children of a cell is or will be @@ -12278,13 +12260,13 @@ bool Triangulation::prepare_coarsening_and_refinement () // eliminate_unrefined_islands, eliminate_refined_inner_islands // and eliminate_refined_boundary_islands will be fulfilled // automatically and do not need to be enforced separately. - // + // // 6/ take care of the requirement that no double refinement is done // at each face - // + // // 7/ take care that no double refinement is done at each line in 3d // or higher dimensions. - // + // // 8/ make sure that all children of each cell are either flagged // for coarsening or none of the children is // @@ -12448,7 +12430,7 @@ bool Triangulation::prepare_coarsening_and_refinement () } ////////////////////////////////////// - // STEP 3: + // STEP 3: // limit the level difference of neighboring cells at each // vertex. // @@ -13290,21 +13272,24 @@ Triangulation::remove_refinement_listener (RefinementListener &li } template <> -const Manifold<2,1> & Triangulation<2, 1>::get_manifold(const types::manifold_id) const { +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 { +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 { +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; diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 69a5ed21a2..6cea86cb86 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1110,13 +1110,13 @@ set_all_manifold_ids (const types::manifold_id manifold_ind) const for (unsigned int c=0; cn_children(); ++c) this->child(c)->set_all_manifold_ids (manifold_ind); - // for hexes also set manifold_id - // of bounding quads and lines - - // Six bonding quads - for(unsigned int i=0; i<6; ++i) + // for hexes also set manifold_id + // of bounding quads and lines + + // Six bonding quads + for (unsigned int i=0; i<6; ++i) this->quad(i)->set_manifold_id(manifold_ind); - // Twelve bounding lines + // Twelve bounding lines for (unsigned int i=0; i<12; ++i) this->line(i)->set_manifold_id (manifold_ind); } diff --git a/source/grid/tria_objects.cc b/source/grid/tria_objects.cc index a61c028abb..3cbe492f1a 100644 --- a/source/grid/tria_objects.cc +++ b/source/grid/tria_objects.cc @@ -202,8 +202,8 @@ namespace internal manifold_id.reserve (new_size); manifold_id.insert (manifold_id.end(), - new_size-manifold_id.size(), - numbers::flat_manifold_id); + new_size-manifold_id.size(), + numbers::flat_manifold_id); } @@ -282,8 +282,8 @@ namespace internal manifold_id.reserve (new_size); manifold_id.insert (manifold_id.end(), - new_size-manifold_id.size(), - numbers::flat_manifold_id); + new_size-manifold_id.size(), + numbers::flat_manifold_id); user_data.reserve (new_size); user_data.resize (new_size); @@ -513,7 +513,7 @@ namespace internal MemoryConsumption::memory_consumption (used) + MemoryConsumption::memory_consumption (user_flags) + MemoryConsumption::memory_consumption (boundary_or_material_id) + - MemoryConsumption::memory_consumption (manifold_id) + + MemoryConsumption::memory_consumption (manifold_id) + MemoryConsumption::memory_consumption (refinement_cases) + user_data.capacity() * sizeof(UserData) + sizeof(user_data)); }