From d8f5f5d847aecc87367cd480b8740b2edf270d51 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 23 Aug 2013 15:48:29 +0000 Subject: [PATCH] Change MappingQ support points from equispaced points into Gauss-Lobatto points. Let's see next if there are more tests failing outside fe. git-svn-id: https://svn.dealii.org/trunk@30450 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/fe/mapping_c1.h | 67 +- deal.II/include/deal.II/fe/mapping_q.h | 400 ++++------- .../include/deal.II/fe/mapping_q_eulerian.h | 88 +-- deal.II/include/deal.II/grid/tria.h | 2 +- deal.II/include/deal.II/grid/tria_boundary.h | 509 +++++--------- deal.II/source/fe/mapping_c1.cc | 114 +--- deal.II/source/fe/mapping_q.cc | 635 +++++------------- deal.II/source/fe/mapping_q_eulerian.cc | 78 +-- deal.II/source/grid/tria_boundary.cc | 92 ++- deal.II/source/grid/tria_boundary_lib.cc | 226 +++---- tests/fe/jacobians/cmp/generic | 40 +- tests/fe/mapping/cmp/generic | 68 +- tests/fe/mapping_q_eulerian/cmp/generic | 2 +- tests/fe/mapping_real_to_unit_q4_sphere.cc | 3 + .../cmp/generic | 1 - 15 files changed, 788 insertions(+), 1537 deletions(-) diff --git a/deal.II/include/deal.II/fe/mapping_c1.h b/deal.II/include/deal.II/fe/mapping_c1.h index 85f900bf28..62bbafecda 100644 --- a/deal.II/include/deal.II/fe/mapping_c1.h +++ b/deal.II/include/deal.II/fe/mapping_c1.h @@ -37,7 +37,7 @@ DEAL_II_NAMESPACE_OPEN * * To use this class, make sure that the * Boundary::@p get_normals_at_vertices function is implemented - * for the users boundary object. + * for the user's boundary object. * * For more information about the spacedim template parameter * check the documentation of FiniteElement or the one of @@ -50,71 +50,46 @@ class MappingC1 : public MappingQ { public: /** - * Constructor. Pass the fixed - * degree @p 3 down to the base - * class, as a cubic mapping - * suffices to generate a - * continuous mapping of the - * boundary. + * Constructor. Pass the fixed degree @p 3 down to the base class, as a + * cubic mapping suffices to generate a continuous mapping of the boundary. */ MappingC1 (); /** - * Return a pointer to a copy of the - * present object. The caller of this - * copy then assumes ownership of it. + * Return a pointer to a copy of the present object. The caller of this copy + * then assumes ownership of it. */ virtual Mapping *clone () const; protected: /** - * For dim=2,3. Append the - * support points of all shape - * functions located on bounding - * lines to the vector - * @p a. Points located on the - * line but on vertices are not - * included. + * For dim=2,3. Append the support points of all shape functions + * located on bounding lines to the vector @p a. Points located on the line + * but on vertices are not included. * - * Needed by the - * compute_support_points_simple(laplace) - * functions. For dim=1 this - * function is empty. + * Needed by the compute_support_points_simple(laplace) + * functions. For dim=1 this function is empty. * - * This function chooses the - * respective points not such - * that they are interpolating - * the boundary (as does the base - * class), but rather such that - * the resulting cubic mapping is - * a continuous one. + * This function chooses the respective points not such that they are + * interpolating the boundary (as does the base class), but rather such that + * the resulting cubic mapping is a continuous one. */ virtual void add_line_support_points (const typename Triangulation::cell_iterator &cell, std::vector > &a) const; /** - * For dim=3. Append the - * support points of all shape - * functions located on bounding - * faces (quads in 3d) to the - * vector @p a. Points located - * on the line but on vertices - * are not included. + * For dim=3. Append the support points of all shape functions + * located on bounding faces (quads in 3d) to the vector @p a. Points + * located on the line but on vertices are not included. * - * Needed by the - * @p compute_support_points_laplace - * function. For dim=1 and 2 - * this function is empty. + * Needed by the @p compute_support_points_laplace function. For + * dim=1 and 2 this function is empty. * - * This function chooses the - * respective points not such - * that they are interpolating - * the boundary (as does the base - * class), but rather such that - * the resulting cubic mapping is - * a continuous one. + * This function chooses the respective points not such that they are + * interpolating the boundary (as does the base class), but rather such that + * the resulting cubic mapping is a continuous one. */ virtual void add_quad_support_points(const typename Triangulation::cell_iterator &cell, diff --git a/deal.II/include/deal.II/fe/mapping_q.h b/deal.II/include/deal.II/fe/mapping_q.h index 484963727b..28e4aefcbf 100644 --- a/deal.II/include/deal.II/fe/mapping_q.h +++ b/deal.II/include/deal.II/fe/mapping_q.h @@ -32,9 +32,11 @@ template class TensorProductPolynomials; /*@{*/ /** - * Mapping class that uses Qp-mappings on boundary cells. The mapping - * shape functions make use of tensor product polynomials with - * equidistant (on the unit cell) support points. + * Mapping class that uses Qp-mappings on boundary cells. The mapping shape + * functions make use of tensor product polynomials with unit cell support + * points equal to the points of the Gauss-Lobatto quadrature formula. These + * points give a well-conditioned interpolation also for very high orders and + * are therefore preferred over equidistant support points. * * For more details about Qp-mappings, see the `mapping' report at * deal.II/doc/reports/mapping_q/index.html in the `Reports' @@ -44,6 +46,10 @@ template class TensorProductPolynomials; * check the documentation of FiniteElement or the one of * Triangulation. * + * @note Since the boundary description is closely tied to the unit cell + * support points, new boundary descriptions need to explicitly use the + * Gauss-Lobatto points. + * * @author Ralf Hartmann, 2000, 2001, 2005; Guido Kanschat 2000, 2001 */ template @@ -51,41 +57,25 @@ class MappingQ : public MappingQ1 { public: /** - * Constructor. @p p gives the - * degree of mapping polynomials - * on boundary cells. + * Constructor. @p p gives the degree of mapping polynomials on boundary + * cells. * - * The second argument determines - * whether the higher order - * mapping should also be used on - * interior cells. If its value - * is false (the - * default), the a lower-order - * mapping is used in the - * interior. This is sufficient - * for most cases where higher - * order mappings are only used - * to better approximate the - * boundary. In that case, cells - * bounded by straight lines are - * acceptable in the - * interior. However, there are - * cases where one would also - * like to use a higher order - * mapping in the interior. The - * MappingQEulerian class is one - * such case. + * The second argument determines whether the higher order mapping should + * also be used on interior cells. If its value is false (the + * default), the a lower-order mapping is used in the interior. This is + * sufficient for most cases where higher order mappings are only used to + * better approximate the boundary. In that case, cells bounded by straight + * lines are acceptable in the interior. However, there are cases where one + * would also like to use a higher order mapping in the interior. The + * MappingQEulerian class is one such case. */ MappingQ (const unsigned int p, const bool use_mapping_q_on_all_cells = false); /** - * Copy constructor. Performs a - * deep copy, i.e. duplicates - * what #tensor_pols points to - * instead of simply copying the - * #tensor_pols pointer as done - * by a default copy constructor. + * Copy constructor. Performs a deep copy, i.e. duplicates what #tensor_pols + * points to instead of simply copying the #tensor_pols pointer as done by a + * default copy constructor. */ MappingQ (const MappingQ &mapping); @@ -95,10 +85,8 @@ public: virtual ~MappingQ (); /** - * Transforms the point @p p on - * the unit cell to the point - * @p p_real on the real cell - * @p cell and returns @p p_real. + * Transforms the point @p p on the unit cell to the point @p p_real on the + * real cell @p cell and returns @p p_real. */ virtual Point transform_unit_to_real_cell ( @@ -106,46 +94,25 @@ public: const Point &p) const; /** - * Transforms the point @p p on - * the real cell to the point - * @p p_unit on the unit cell - * @p cell and returns @p p_unit. + * Transforms the point @p p on the real cell to the point @p p_unit on the + * unit cell @p cell and returns @p p_unit. * - * Uses Newton iteration and the - * @p transform_unit_to_real_cell - * function. + * Uses Newton iteration and the @p transform_unit_to_real_cell function. * - * In the codimension one case, - * this function returns the - * normal projection of the real - * point @p p on the curve or - * surface identified by the @p - * cell. + * In the codimension one case, this function returns the normal projection + * of the real point @p p on the curve or surface identified by the @p cell. * - * @note Polynomial mappings from - * the reference (unit) cell coordinates - * to the coordinate system of a real - * cell are not always invertible if - * the point for which the inverse - * mapping is to be computed lies - * outside the cell's boundaries. - * In such cases, the current function - * may fail to compute a point on - * the reference cell whose image - * under the mapping equals the given - * point @p p. If this is the case - * then this function throws an - * exception of type - * Mapping::ExcTransformationFailed . - * Whether the given point @p p lies - * outside the cell can therefore be - * determined by checking whether the - * return reference coordinates lie - * inside of outside the reference - * cell (e.g., using - * GeometryInfo::is_inside_unit_cell) - * or whether the exception mentioned - * above has been thrown. + * @note Polynomial mappings from the reference (unit) cell coordinates to + * the coordinate system of a real cell are not always invertible if the + * point for which the inverse mapping is to be computed lies outside the + * cell's boundaries. In such cases, the current function may fail to + * compute a point on the reference cell whose image under the mapping + * equals the given point @p p. If this is the case then this function + * throws an exception of type Mapping::ExcTransformationFailed . Whether + * the given point @p p lies outside the cell can therefore be determined by + * checking whether the return reference coordinates lie inside of outside + * the reference cell (e.g., using GeometryInfo::is_inside_unit_cell) or + * whether the exception mentioned above has been thrown. */ virtual Point transform_real_to_unit_cell ( @@ -172,23 +139,20 @@ public: const MappingType type) const; /** - * Return the degree of the - * mapping, i.e. the value which - * was passed to the constructor. + * Return the degree of the mapping, i.e. the value which was passed to the + * constructor. */ unsigned int get_degree () const; /** - * Return a pointer to a copy of the - * present object. The caller of this - * copy then assumes ownership of it. + * Return a pointer to a copy of the present object. The caller of this copy + * then assumes ownership of it. */ virtual Mapping *clone () const; /** - * Storage for internal data of - * Q_degree transformation. + * Storage for internal data of Q_degree transformation. */ class InternalData : public MappingQ1::InternalData { @@ -200,50 +164,35 @@ public: /** - * Return an estimate (in - * bytes) or the memory - * consumption of this - * object. + * Return an estimate (in bytes) or the memory consumption of this object. */ virtual std::size_t memory_consumption () const; /** - * Unit normal vectors. Used - * for the alternative - * computation of the normal - * vectors. See doc of the - * @p alternative_normals_computation - * flag. + * Unit normal vectors. Used for the alternative computation of the normal + * vectors. See doc of the @p alternative_normals_computation flag. * - * Filled (hardcoded) once in - * @p get_face_data. + * Filled (hardcoded) once in @p get_face_data. */ std::vector > > unit_normals; /** - * Flag that is set by the - * fill_fe_[[sub]face]_values - * function. + * Flag that is set by the fill_fe_[[sub]face]_values function. * - * If this flag is @p true - * we are on an interior cell - * and the - * @p mapping_q1_data is - * used. + * If this flag is @p true we are on an interior cell and the @p + * mapping_q1_data is used. */ bool use_mapping_q1_on_current_cell; /** - * On interior cells - * @p MappingQ1 is used. + * On interior cells @p MappingQ1 is used. */ typename MappingQ1::InternalData mapping_q1_data; }; protected: /** - * Implementation of the interface in - * Mapping. + * Implementation of the interface in Mapping. */ virtual void fill_fe_values (const typename Triangulation::cell_iterator &cell, @@ -258,8 +207,7 @@ protected: CellSimilarity::Similarity &cell_similarity) const ; /** - * Implementation of the interface in - * Mapping. + * Implementation of the interface in Mapping. */ virtual void fill_fe_face_values (const typename Triangulation::cell_iterator &cell, @@ -272,8 +220,7 @@ protected: typename std::vector > &normal_vectors) const ; /** - * Implementation of the interface in - * Mapping. + * Implementation of the interface in Mapping. */ virtual void fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, @@ -287,55 +234,32 @@ protected: typename std::vector > &normal_vectors) const ; /** - * For dim=2,3. Append the - * support points of all shape - * functions located on bounding - * lines to the vector - * @p a. Points located on the - * line but not on vertices are not - * included. + * For dim=2,3. Append the support points of all shape functions + * located on bounding lines to the vector @p a. Points located on the line + * but not on vertices are not included. * - * Needed by the - * @p compute_support_points_laplace - * function . For dim=1 this - * function is empty. + * Needed by the @p compute_support_points_laplace function . For + * dim=1 this function is empty. * - * This function is made virtual - * in order to allow derived - * classes to choose shape - * function support points - * differently than the present - * class, which chooses the - * points as interpolation points - * on the boundary. + * This function is made virtual in order to allow derived classes to choose + * shape function support points differently than the present class, which + * chooses the points as interpolation points on the boundary. */ virtual void add_line_support_points (const typename Triangulation::cell_iterator &cell, std::vector > &a) const; /** - * For dim=3. Append the - * support points of all shape - * functions located on bounding - * faces (quads in 3d) to the - * vector @p a. Points located - * on the quad but not on vertices - * are not included. + * For dim=3. Append the support points of all shape functions + * located on bounding faces (quads in 3d) to the vector @p a. Points + * located on the quad but not on vertices are not included. * - * Needed by the - * @p compute_support_points_laplace - * function. For dim=1 and - * dim=2 this function is - * empty. + * Needed by the @p compute_support_points_laplace function. For + * dim=1 and dim=2 this function is empty. * - * This function is made virtual - * in order to allow derived - * classes to choose shape - * function support points - * differently than the present - * class, which chooses the - * points as interpolation points - * on the boundary. + * This function is made virtual in order to allow derived classes to choose + * shape function support points differently than the present class, which + * chooses the points as interpolation points on the boundary. */ virtual void add_quad_support_points(const typename Triangulation::cell_iterator &cell, @@ -359,145 +283,99 @@ private: const Quadrature& quadrature) const; /** - * Compute shape values and/or - * derivatives. + * Compute shape values and/or derivatives. */ virtual void compute_shapes_virtual (const std::vector > &unit_points, typename MappingQ1::InternalData &data) const; /** - * This function is needed by the - * constructor of MappingQ - * for dim= 2 and 3. + * This function is needed by the constructor of + * MappingQ for dim= 2 and 3. * - * For degree<4 this function - * sets the - * @p laplace_on_quad_vector to - * the hardcoded data. For - * degree>=4 and MappingQ<2> - * this vector is computed. + * For degree<4 this function sets the @p laplace_on_quad_vector to + * the hardcoded data. For degree>=4 and MappingQ<2> this vector is + * computed. * - * For the definition of the - * @p laplace_on_quad_vector - * please refer to equation (8) - * of the `mapping' report. + * For the definition of the @p laplace_on_quad_vector please refer to + * equation (8) of the `mapping' report. */ void set_laplace_on_quad_vector(Table<2,double> &loqvs) const; /** - * This function is needed by the - * constructor of MappingQ<3>. + * This function is needed by the constructor of MappingQ<3>. * - * For degree==2 this function - * sets the - * @p laplace_on_hex_vector to - * the hardcoded data. For - * degree>2 this vector is - * computed. + * For degree==2 this function sets the @p laplace_on_hex_vector to + * the hardcoded data. For degree>2 this vector is computed. * - * For the definition of the - * @p laplace_on_hex_vector - * please refer to equation (8) - * of the `mapping' report. + * For the definition of the @p laplace_on_hex_vector please refer to + * equation (8) of the `mapping' report. */ void set_laplace_on_hex_vector(Table<2,double> &lohvs) const; /** - * Computes the - * laplace_on_quad(hex)_vector. + * Computes the laplace_on_quad(hex)_vector. * - * Called by the - * set_laplace_on_quad(hex)_vector - * functions if the data is not - * yet hardcoded. + * Called by the set_laplace_on_quad(hex)_vector functions if the + * data is not yet hardcoded. * - * For the definition of the - * laplace_on_quad(hex)_vector - * please refer to equation (8) - * of the `mapping' report. + * For the definition of the laplace_on_quad(hex)_vector please + * refer to equation (8) of the `mapping' report. */ void compute_laplace_vector(Table<2,double> &lvs) const; /** - * Takes a - * laplace_on_hex(quad)_vector - * and applies it to the vector - * @p a to compute the inner - * support points as a linear - * combination of the exterior - * points. + * Takes a laplace_on_hex(quad)_vector and applies it to the vector + * @p a to compute the inner support points as a linear combination of the + * exterior points. * - * The vector @p a initially - * contains the locations of the - * @p n_outer points, the - * @p n_inner computed inner - * points are appended. + * The vector @p a initially contains the locations of the @p n_outer + * points, the @p n_inner computed inner points are appended. * - * See equation (7) of the - * `mapping' report. + * See equation (7) of the `mapping' report. */ void apply_laplace_vector(const Table<2,double> &lvs, std::vector > &a) const; /** - * Computes the support points of - * the mapping. + * Computes the support points of the mapping. */ virtual void compute_mapping_support_points( const typename Triangulation::cell_iterator &cell, std::vector > &a) const; /** - * Computes all support points of - * the mapping shape - * functions. The inner support - * points (ie. support points in - * quads for 2d, in hexes for 3d) - * are computed using the - * solution of a Laplace equation - * with the position of the outer - * support points as boundary - * values, in order to make the - * transformation as smooth as - * possible. + * Computes all support points of the mapping shape functions. The inner + * support points (ie. support points in quads for 2d, in hexes for 3d) are + * computed using the solution of a Laplace equation with the position of + * the outer support points as boundary values, in order to make the + * transformation as smooth as possible. */ void compute_support_points_laplace( const typename Triangulation::cell_iterator &cell, std::vector > &a) const; /** - * Needed by the - * @p laplace_on_quad function - * (for dim==2). Filled by the - * constructor. + * Needed by the @p laplace_on_quad function (for dim==2). Filled + * by the constructor. * * Sizes: - * laplace_on_quad_vector.size()= - * number of inner - * unit_support_points - * laplace_on_quad_vector[i].size()= - * number of outer - * unit_support_points, i.e. - * unit_support_points on the - * boundary of the quad + * laplace_on_quad_vector.size()= number of inner unit_support_points + * laplace_on_quad_vector[i].size()= number of outer unit_support_points, + * i.e. unit_support_points on the boundary of the quad * - * For the definition of this - * vector see equation (8) of the - * `mapping' report. + * For the definition of this vector see equation (8) of the `mapping' + * report. */ Table<2,double> laplace_on_quad_vector; /** - * Needed by the - * @p laplace_on_hex function - * (for dim==3). Filled by the - * constructor. + * Needed by the @p laplace_on_hex function (for dim==3). Filled by + * the constructor. * - * For the definition of this - * vector see equation (8) of the - * `mapping' report. + * For the definition of this vector see equation (8) of the `mapping' + * report. */ Table<2,double> laplace_on_hex_vector; @@ -509,69 +387,57 @@ private: << "laplace_vector not set for degree=" << arg1 << "."); /** - * Degree @p p of the - * polynomials used as shape - * functions for the Qp mapping + * Degree @p p of the polynomials used as shape functions for the Qp mapping * of cells at the boundary. */ const unsigned int degree; /** - * Number of inner mapping shape - * functions. + * Number of inner mapping shape functions. */ const unsigned int n_inner; /** - * Number of mapping shape - * functions on the boundary. + * Number of mapping shape functions on the boundary. */ const unsigned int n_outer; /** - * Pointer to the - * @p dim-dimensional tensor - * product polynomials used as - * shape functions for the Qp - * mapping of cells at the - * boundary. + * Pointer to the @p dim-dimensional tensor product polynomials used as + * shape functions for the Qp mapping of cells at the boundary. */ const TensorProductPolynomials *tensor_pols; /** - * Number of the Qp tensor - * product shape functions. + * Number of the Qp tensor product shape functions. */ const unsigned int n_shape_functions; /** - * Mapping from lexicographic to - * to the Qp shape function - * numbering. Its size is - * @p dofs_per_cell. + * Mapping from lexicographic to to the Qp shape function numbering. Its + * size is @p dofs_per_cell. */ const std::vector renumber; /** - * If this flag is set @p true - * then @p MappingQ is used on - * all cells, not only on - * boundary cells. + * If this flag is set @p true then @p MappingQ is used on all cells, not + * only on boundary cells. */ const bool use_mapping_q_on_all_cells; /** - * An FE_Q object which is only needed in - * 3D, since it knows how to reorder shape - * functions/DoFs on non-standard - * faces. This is used to reorder support - * points in the same way. We could make - * this a pointer to prevent construction - * in 1D and 2D, but since memory and time - * requirements are not particularly high - * this seems unnecessary at the moment. + * An FE_Q object which is only needed in 3D, since it knows how to reorder + * shape functions/DoFs on non-standard faces. This is used to reorder + * support points in the same way. We could make this a pointer to prevent + * construction in 1D and 2D, but since memory and time requirements are not + * particularly high this seems unnecessary at the moment. */ const FE_Q feq; + + /** + * Declare other MappingQ classes friends. + */ + template friend class MappingQ; }; /*@}*/ diff --git a/deal.II/include/deal.II/fe/mapping_q_eulerian.h b/deal.II/include/deal.II/fe/mapping_q_eulerian.h index 4d4ba9cc9c..f0a79f1066 100644 --- a/deal.II/include/deal.II/fe/mapping_q_eulerian.h +++ b/deal.II/include/deal.II/fe/mapping_q_eulerian.h @@ -95,25 +95,14 @@ class MappingQEulerian : public MappingQ { public: /** - * Constructor. The first argument is - * the polynomical degree of the desired - * Qp mapping. It then takes a - * Vector & to specify the - * transformation of the domain - * from the reference to - * the current configuration. - * The organization of the - * elements in the @p Vector - * must follow the concept how - * deal.II stores solutions that - * are associated to a - * triangulation. This is - * automatically the case if the - * @p Vector represents the - * solution of the previous step - * of a nonlinear problem. - * Alternatively, the @p Vector - * can be initialized by + * Constructor. The first argument is the polynomical degree of the desired + * Qp mapping. It then takes a Vector & to specify the + * transformation of the domain from the reference to the current + * configuration. The organization of the elements in the @p Vector must + * follow the concept how deal.II stores solutions that are associated to a + * triangulation. This is automatically the case if the @p Vector + * represents the solution of the previous step of a nonlinear problem. + * Alternatively, the @p Vector can be initialized by * DoFAccessor::set_dof_values(). */ @@ -122,20 +111,16 @@ public: const DoFHandler &euler_dof_handler); /** - * Return a pointer to a copy of the - * present object. The caller of this - * copy then assumes ownership of it. + * Return a pointer to a copy of the present object. The caller of this copy + * then assumes ownership of it. */ virtual Mapping *clone () const; /** - * Always returns @p false because - * MappingQ1Eulerian does not in general - * preserve vertex locations (unless the - * translation vector happens to provide - * for zero displacements at vertex - * locations). + * Always returns @p false because MappingQ1Eulerian does not in general + * preserve vertex locations (unless the translation vector happens to + * provide for zero displacements at vertex locations). */ bool preserves_vertex_locations () const; @@ -146,10 +131,8 @@ public: protected: /** - * Implementation of the interface in - * MappingQ. Overrides the function in - * the base class, since we cannot use - * any cell similarity for this class. + * Implementation of the interface in MappingQ. Overrides the function in + * the base class, since we cannot use any cell similarity for this class. */ virtual void fill_fe_values (const typename Triangulation::cell_iterator &cell, @@ -164,36 +147,29 @@ protected: CellSimilarity::Similarity &cell_similarity) const; /** - * Reference to the vector of - * shifts. + * Reference to the vector of shifts. */ SmartPointer > euler_vector; /** - * Pointer to the DoFHandler to - * which the mapping vector is - * associated. + * Pointer to the DoFHandler to which the mapping vector is associated. */ - SmartPointer,MappingQEulerian > euler_dof_handler; private: /** - * Special quadrature rule used - * to define the support points - * in the reference configuration. + * Special quadrature rule used to define the support points in the + * reference configuration. */ class SupportQuadrature : public Quadrature { public: /** - * Constructor, with an argument - * defining the desired polynomial - * degree. + * Constructor, with an argument defining the desired polynomial degree. */ SupportQuadrature (const unsigned int map_degree); @@ -201,36 +177,26 @@ private: }; /** - * A member variable holding the - * quadrature points in the right - * order. + * A member variable holding the quadrature points in the right order. */ const SupportQuadrature support_quadrature; /** - * FEValues object used to query the - * the given finite element field - * at the support points in the - * reference configuration. + * FEValues object used to query the the given finite element field at the + * support points in the reference configuration. * - * The variable is marked as - * mutable since we have to call - * FEValues::reinit from - * compute_mapping_support_points, - * a function that is 'const'. + * The variable is marked as mutable since we have to call FEValues::reinit + * from compute_mapping_support_points, a function that is 'const'. */ mutable FEValues fe_values; /** - * A variable to guard access to - * the fe_values variable. + * A variable to guard access to the fe_values variable. */ mutable Threads::Mutex fe_values_mutex; /** - * Compute the positions of the - * support points in the current - * configuration + * Compute the positions of the support points in the current configuration */ virtual void compute_mapping_support_points( const typename Triangulation::cell_iterator &cell, diff --git a/deal.II/include/deal.II/grid/tria.h b/deal.II/include/deal.II/grid/tria.h index 72b71db88f..b478da2009 100644 --- a/deal.II/include/deal.II/grid/tria.h +++ b/deal.II/include/deal.II/grid/tria.h @@ -1566,7 +1566,7 @@ public: virtual void create_notification (const Triangulation &tria); - } DEAL_II_DEPRECATED; + }; /** * A structure that is used as an diff --git a/deal.II/include/deal.II/grid/tria_boundary.h b/deal.II/include/deal.II/grid/tria_boundary.h index ba83065200..ca948550fc 100644 --- a/deal.II/include/deal.II/grid/tria_boundary.h +++ b/deal.II/include/deal.II/grid/tria_boundary.h @@ -22,6 +22,8 @@ #include #include +#include +#include #include #include @@ -85,120 +87,75 @@ class Boundary : public Subscriptor public: /** - * Type keeping information about - * the normals at the vertices of - * a face of a cell. Thus, there - * are - * GeometryInfo::vertices_per_face - * normal vectors, that define - * the tangent spaces of the - * boundary at the vertices. Note - * that the vectors stored in - * this object are not required - * to be normalized, nor to - * actually point outward, as one - * often will only want to check - * for orthogonality to define - * the tangent plane; if a - * function requires the normals - * to be normalized, then it must - * do so itself. + * Type keeping information about the normals at the vertices of a face of a + * cell. Thus, there are GeometryInfo::vertices_per_face + * normal vectors, that define the tangent spaces of the boundary at the + * vertices. Note that the vectors stored in this object are not required to + * be normalized, nor to actually point outward, as one often will only want + * to check for orthogonality to define the tangent plane; if a function + * requires the normals to be normalized, then it must do so itself. * - * For obvious reasons, this - * type is not useful in 1d. + * For obvious reasons, this type is not useful in 1d. */ typedef Tensor<1,spacedim> FaceVertexNormals[GeometryInfo::vertices_per_face]; /** - * Destructor. Does nothing here, but - * needs to be declared to make it + * Destructor. Does nothing here, but needs to be declared to make it * virtual. */ virtual ~Boundary (); /** - * Return the point which shall - * become the new middle vertex - * of the two children of a - * regular line. In 2D, this line - * is a line at the boundary, - * while in 3d, it is bounding a - * face at the boundary (the - * lines therefore is also on the - * boundary). + * Return the point which shall become the new middle vertex of the two + * children of a regular line. In 2D, this line is a line at the boundary, + * while in 3d, it is bounding a face at the boundary (the lines therefore + * is also on the boundary). */ virtual Point get_new_point_on_line (const typename Triangulation::line_iterator &line) const = 0; /** - * Return the point which shall - * become the common point of the - * four children of a quad at the - * boundary in three or more - * spatial dimensions. This - * function therefore is only - * useful in at least three - * dimensions and should not be - * called for lower dimensions. + * Return the point which shall become the common point of the four children + * of a quad at the boundary in three or more spatial dimensions. This + * function therefore is only useful in at least three dimensions and should + * not be called for lower dimensions. * - * This function is called after - * the four lines bounding the - * given @p quad are refined, so - * you may want to use the - * information provided by - * quad->line(i)->child(j), - * i=0...3, j=0,1. + * This function is called after the four lines bounding the given @p quad + * are refined, so you may want to use the information provided by + * quad->line(i)->child(j), i=0...3, j=0,1. * - * Because in 2D, this function - * is not needed, it is not made - * pure virtual, to avoid the - * need to overload it. The - * default implementation throws - * an error in any case, however. + * Because in 2D, this function is not needed, it is not made pure virtual, + * to avoid the need to overload it. The default implementation throws an + * error in any case, however. */ virtual Point get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; /** - * Depending on dim=2 or - * dim=3 this function - * calls the - * get_new_point_on_line or the - * get_new_point_on_quad - * function. It throws an - * exception for - * dim=1. This wrapper - * allows dimension independent + * Depending on dim=2 or dim=3 this function calls the + * get_new_point_on_line or the get_new_point_on_quad function. It throws an + * exception for dim=1. This wrapper allows dimension independent * programming. */ Point get_new_point_on_face (const typename Triangulation::face_iterator &face) const; /** - * Return equally spaced - * intermediate points on a line. + * Return intermediate points on a line spaced according to the interior + * support points of the 1D Gauss-Lobatto quadrature formula. * - * The number of points requested - * is given by the size of the - * vector @p points. It is the - * task of the derived classes to - * arrange the points in + * The number of points requested is given by the size of the vector @p + * points. It is the task of the derived classes to arrange the points in * approximately equal distances. * - * This function is called by the - * @p MappingQ class. This - * happens on each face line of a - * cells that has got at least - * one boundary line. + * This function is called by the @p MappingQ class. This happens on each + * face line of a cells that has got at least one boundary line. * - * As this function is not needed - * for @p MappingQ1, it is not - * made pure virtual, to avoid - * the need to overload it. The - * default implementation throws - * an error in any case, however. + * As this function is not needed for @p MappingQ1, it is not made pure + * virtual, to avoid the need to overload it. The default implementation + * throws an error in any case, however. */ virtual void @@ -206,34 +163,23 @@ public: std::vector > &points) const; /** - * Return equally spaced - * intermediate points on a - * boundary quad. + * Return intermediate points on a line spaced according to the tensor + * product of the interior support points of the 1D Gauss-Lobatto quadrature + * formula. * - * The number of points requested - * is given by the size of the - * vector @p points. It is - * required that this number is a - * square of another integer, - * i.e. n=points.size()=m*m. It - * is the task of the derived - * classes to arrange the points - * such they split the quad into - * (m+1)(m+1) approximately - * equal-sized subquads. + * The number of points requested is given by the size of the vector @p + * points. It is required that this number is a square of another integer, + * i.e. n=points.size()=m*m. It is the task of the derived classes + * to arrange the points such they split the quad into (m+1)(m+1) + * approximately equal-sized subquads. * - * This function is called by the - * MappingQ<3> class. This - * happens each face quad of - * cells in 3d that has got at - * least one boundary face quad. + * This function is called by the MappingQ<3> class. This happens + * each face quad of cells in 3d that has got at least one boundary face + * quad. * - * As this function is not needed - * for @p MappingQ1, it is not - * made pure virtual, to avoid - * the need to overload it. The - * default implementation throws - * an error in any case, however. + * As this function is not needed for @p MappingQ1, it is not made pure + * virtual, to avoid the need to overload it. The default implementation + * throws an error in any case, however. */ virtual void @@ -241,56 +187,36 @@ public: std::vector > &points) const; /** - * Depending on dim=2 or - * dim=3 this function - * calls the - * get_intermediate_points_on_line - * or the - * get_intermediate_points_on_quad - * function. It throws an - * exception for - * dim=1. This wrapper - * allows dimension independent - * programming. + * Depending on dim=2 or dim=3 this function calls the + * get_intermediate_points_on_line or the get_intermediate_points_on_quad + * function. It throws an exception for dim=1. This wrapper allows + * dimension independent programming. */ void get_intermediate_points_on_face (const typename Triangulation::face_iterator &face, std::vector > &points) const; /** - * Return the normal vector to the surface - * at the point p. If p is not in fact - * on the surface, but only closeby, - * try to return something reasonable, - * for example the normal vector - * at the surface point closest to p. - * (The point p will in fact not normally - * lie on the actual surface, but rather - * be a quadrature point mapped by some - * polynomial mapping; the mapped surface, - * however, will not usually coincide with - * the actual surface.) + * Return the normal vector to the surface at the point p. If p is not in + * fact on the surface, but only close-by, try to return something + * reasonable, for example the normal vector at the surface point closest to + * p. (The point p will in fact not normally lie on the actual surface, but + * rather be a quadrature point mapped by some polynomial mapping; the + * mapped surface, however, will not usually coincide with the actual + * surface.) * - * The face iterator gives an indication - * which face this function is supposed - * to compute the normal vector for. - * This is useful if the boundary of - * the domain is composed of different - * nondifferential pieces (for example - * when using the StraightBoundary class - * to approximate a geometry that is - * completely described by the coarse mesh, - * with piecewise (bi-)linear components - * between the vertices, but where the - * boundary may have a kink at the vertices - * itself). + * The face iterator gives an indication which face this function is + * supposed to compute the normal vector for. This is useful if the + * boundary of the domain is composed of different nondifferential pieces + * (for example when using the StraightBoundary class to approximate a + * geometry that is completely described by the coarse mesh, with piecewise + * (bi-)linear components between the vertices, but where the boundary may + * have a kink at the vertices itself). * - * @note Implementations of this function - * should be able to assume that the point p - * lies within or close to the face described by the - * first argument. In turn, callers of this - * function should ensure that this is - * in fact the case. + * @note Implementations of this function should be able to assume that the + * point p lies within or close to the face described by the first + * argument. In turn, callers of this function should ensure that this is in + * fact the case. */ virtual Tensor<1,spacedim> @@ -298,32 +224,17 @@ public: const Point &p) const; /** - * Compute the normal vectors to - * the boundary at each vertex of - * the given face. It is not - * required that the normal - * vectors be normed - * somehow. Neither is it - * required that the normals - * actually point outward. + * Compute the normal vectors to the boundary at each vertex of the given + * face. It is not required that the normal vectors be normed + * somehow. Neither is it required that the normals actually point outward. * - * This function is - * needed to compute data for C1 - * mappings. The default - * implementation is to throw an - * error, so you need not - * overload this function in case - * you do not intend to use C1 - * mappings. + * This function is needed to compute data for C1 mappings. The default + * implementation is to throw an error, so you need not overload this + * function in case you do not intend to use C1 mappings. * - * Note that when computing - * normal vectors at a vertex - * where the boundary is not - * differentiable, you have to - * make sure that you compute the - * one-sided limits, i.e. limit - * with respect to points inside - * the given face. + * Note that when computing normal vectors at a vertex where the boundary is + * not differentiable, you have to make sure that you compute the one-sided + * limits, i.e. limit with respect to points inside the given face. */ virtual void @@ -331,35 +242,20 @@ public: FaceVertexNormals &face_vertex_normals) const; /** - * Given a candidate point and a - * line segment characterized by - * the iterator, return a point - * that lies on the surface - * described by this object. This - * function is used in some mesh - * smoothing algorithms that try - * to move around points in order - * to improve the mesh quality - * but need to ensure that points - * that were on the boundary - * remain on the boundary. + * Given a candidate point and a line segment characterized by the iterator, + * return a point that lies on the surface described by this object. This + * function is used in some mesh smoothing algorithms that try to move + * around points in order to improve the mesh quality but need to ensure + * that points that were on the boundary remain on the boundary. * - * If spacedim==1, then the line - * represented by the line - * iterator is the entire space - * (i.e. it is a cell, not a part - * of the boundary), and the - * returned point equals the - * given input point. + * If spacedim==1, then the line represented by the line iterator is the + * entire space (i.e. it is a cell, not a part of the boundary), and the + * returned point equals the given input point. * - * Derived classes do not need to - * implement this function unless - * mesh smoothing algorithms are - * used with a particular - * boundary object. The default - * implementation of this - * function throws an exception - * of type ExcPureFunctionCalled. + * Derived classes do not need to implement this function unless mesh + * smoothing algorithms are used with a particular boundary object. The + * default implementation of this function throws an exception of type + * ExcPureFunctionCalled. */ virtual Point @@ -367,19 +263,12 @@ public: const Point &candidate) const; /** - * Same function as above but for - * a point that is to be - * projected onto the area - * characterized by the given - * quad. + * Same function as above but for a point that is to be projected onto the + * area characterized by the given quad. * - * If spacedim<=2, then the surface - * represented by the quad - * iterator is the entire space - * (i.e. it is a cell, not a part - * of the boundary), and the - * returned point equals the - * given input point. + * If spacedim<=2, then the surface represented by the quad iterator is the + * entire space (i.e. it is a cell, not a part of the boundary), and the + * returned point equals the given input point. */ virtual Point @@ -387,24 +276,40 @@ public: const Point &candidate) const; /** - * Same function as above but for - * a point that is to be - * projected onto the area - * characterized by the given - * quad. + * Same function as above but for a point that is to be projected onto the + * area characterized by the given quad. * - * If spacedim<=3, then the manifold - * represented by the hex - * iterator is the entire space - * (i.e. it is a cell, not a part - * of the boundary), and the - * returned point equals the - * given input point. + * If spacedim<=3, then the manifold represented by the hex iterator is the + * entire space (i.e. it is a cell, not a part of the boundary), and the + * returned point equals the given input point. */ virtual Point project_to_surface (const typename Triangulation::hex_iterator &hex, const Point &candidate) const; + +protected: + /** + * Returns the support points of the Gauss-Lobatto quadrature formula used + * for intermediate points. + * + * @note Since the boundary description is closely tied to the unit cell + * support points of MappingQ, new boundary descriptions need to explicitly + * use these Gauss-Lobatto points and not equidistant points. + */ + const std::vector > & + get_line_support_points (const unsigned int n_intermediate_points) const; + +private: + /** + * Point generator for the intermediate points on a boundary. + */ + mutable std::vector > > points; + + /** + * Mutex for protecting the points array. + */ + mutable Threads::Mutex mutex; }; @@ -429,56 +334,37 @@ class StraightBoundary : public Boundary { public: /** - * Default constructor. Some - * compilers require this for - * some reasons. + * Default constructor. Some compilers require this for some reasons. */ StraightBoundary (); /** - * Let the new point be the - * arithmetic mean of the two - * vertices of the line. + * Let the new point be the arithmetic mean of the two vertices of the line. * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class for more - * information. + * Refer to the general documentation of this class and the documentation of + * the base class for more information. */ virtual Point get_new_point_on_line (const typename Triangulation::line_iterator &line) const; /** - * Let the new point be the - * arithmetic mean of the four - * vertices of this quad and the - * four midpoints of the lines, - * which are already created at - * the time of calling this - * function. + * Let the new point be the arithmetic mean of the four vertices of this + * quad and the four midpoints of the lines, which are already created at + * the time of calling this function. * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class for more - * information. + * Refer to the general documentation of this class and the documentation of + * the base class for more information. */ virtual Point get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; /** - * Gives n=points.size() - * points that splits the - * StraightBoundary line into - * $n+1$ partitions of equal - * lengths. + * Gives n=points.size() points that splits the StraightBoundary + * line into $n+1$ partitions of equal lengths. * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. + * Refer to the general documentation of this class and the documentation of + * the base class. */ virtual void @@ -486,16 +372,11 @@ public: std::vector > &points) const; /** - * Gives n=points.size()=m*m - * points that splits the - * StraightBoundary quad into - * $(m+1)(m+1)$ subquads of equal - * size. + * Gives n=points.size()=m*m points that splits the + * StraightBoundary quad into $(m+1)(m+1)$ subquads of equal size. * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. + * Refer to the general documentation of this class and the documentation of + * the base class. */ virtual void @@ -503,13 +384,10 @@ public: std::vector > &points) const; /** - * Implementation of the function - * declared in the base class. + * Implementation of the function declared in the base class. * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. + * Refer to the general documentation of this class and the documentation of + * the base class. */ virtual Tensor<1,spacedim> @@ -517,14 +395,10 @@ public: const Point &p) const; /** - * Compute the normals to the - * boundary at the vertices of - * the given face. + * Compute the normals to the boundary at the vertices of the given face. * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. + * Refer to the general documentation of this class and the documentation of + * the base class. */ virtual void @@ -532,32 +406,18 @@ public: typename Boundary::FaceVertexNormals &face_vertex_normals) const; /** - * Given a candidate point and a - * line segment characterized by - * the iterator, return a point - * that lies on the surface - * described by this object. This - * function is used in some mesh - * smoothing algorithms that try - * to move around points in order - * to improve the mesh quality - * but need to ensure that points - * that were on the boundary - * remain on the boundary. + * Given a candidate point and a line segment characterized by the iterator, + * return a point that lies on the surface described by this object. This + * function is used in some mesh smoothing algorithms that try to move + * around points in order to improve the mesh quality but need to ensure + * that points that were on the boundary remain on the boundary. * - * The point returned is the - * projection of the candidate - * point onto the line through - * the two vertices of the given - * line iterator. + * The point returned is the projection of the candidate point onto the line + * through the two vertices of the given line iterator. * - * If spacedim==1, then the line - * represented by the line - * iterator is the entire space - * (i.e. it is a cell, not a part - * of the boundary), and the - * returned point equals the - * given input point. + * If spacedim==1, then the line represented by the line iterator is the + * entire space (i.e. it is a cell, not a part of the boundary), and the + * returned point equals the given input point. */ virtual Point @@ -565,26 +425,15 @@ public: const Point &candidate) const; /** - * Same function as above but for - * a point that is to be - * projected onto the area - * characterized by the given - * quad. + * Same function as above but for a point that is to be projected onto the + * area characterized by the given quad. * - * The point returned is the - * projection of the candidate - * point onto the bilinear - * surface spanned by the four - * vertices of the given quad - * iterator. + * The point returned is the projection of the candidate point onto the + * bilinear surface spanned by the four vertices of the given quad iterator. * - * If spacedim<=2, then the surface - * represented by the quad - * iterator is the entire space - * (i.e. it is a cell, not a part - * of the boundary), and the - * returned point equals the - * given input point. + * If spacedim<=2, then the surface represented by the quad iterator is the + * entire space (i.e. it is a cell, not a part of the boundary), and the + * returned point equals the given input point. */ virtual Point @@ -592,26 +441,16 @@ public: const Point &candidate) const; /** - * Same function as above but for - * a point that is to be - * projected onto the area - * characterized by the given - * quad. + * Same function as above but for a point that is to be projected onto the + * area characterized by the given quad. * - * The point returned is the - * projection of the candidate - * point onto the trilinear - * manifold spanned by the eight - * vertices of the given hex + * The point returned is the projection of the candidate point onto the + * trilinear manifold spanned by the eight vertices of the given hex * iterator. * - * If spacedim<=3, then the manifold - * represented by the hex - * iterator is the entire space - * (i.e. it is a cell, not a part - * of the boundary), and the - * returned point equals the - * given input point. + * If spacedim<=3, then the manifold represented by the hex iterator is the + * entire space (i.e. it is a cell, not a part of the boundary), and the + * returned point equals the given input point. */ virtual Point diff --git a/deal.II/source/fe/mapping_c1.cc b/deal.II/source/fe/mapping_c1.cc index 1873d9c737..5a6115220a 100644 --- a/deal.II/source/fe/mapping_c1.cc +++ b/deal.II/source/fe/mapping_c1.cc @@ -53,88 +53,46 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce const unsigned int dim = 2; std::vector > line_points (2); - // loop over each of the lines, - // and if it is at the - // boundary, then first get the - // boundary description and - // second compute the points on - // it. if not at the boundary, - // get the respective points - // from another function + // loop over each of the lines, and if it is at the boundary, then first get + // the boundary description and second compute the points on it. if not at + // the boundary, get the respective points from another function for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) { const Triangulation::line_iterator line = cell->line(line_no); if (line->at_boundary()) { - // first get the normal - // vectors at the two - // vertices of this line - // from the boundary - // description + // first get the normal vectors at the two vertices of this line + // from the boundary description const Boundary &boundary = line->get_triangulation().get_boundary(line->boundary_indicator()); Boundary::FaceVertexNormals face_vertex_normals; boundary.get_normals_at_vertices (line, face_vertex_normals); - // then transform them into - // interpolation points for - // a cubic polynomial + // then transform them into interpolation points for a cubic + // polynomial // - // for this, note that if - // we describe the boundary - // curve as a polynomial in - // tangential coordinate - // @p{t=0..1} (along the - // line) and @p{s} in - // normal direction, then - // the cubic mapping is - // such that @p{s = a*t**3 - // + b*t**2 + c*t + d}, and - // we want to determine the - // interpolation points at - // @p{t=1/3} and - // @p{t=2/3}. Since at - // @p{t=0,1} we want a - // vertex which is actually - // at the boundary, we know - // that @p{d=0} and - // @p{a=-b-c}. As - // side-conditions, we want - // that the derivatives at - // @p{t=0} and @p{t=1}, - // i.e. at the vertices - // match those returned by - // the boundary. We then - // have that - // @p{s(1/3)=1/27(2b+8c)} - // and - // @p{s(2/3)=4/27b+10/27c}. + // for this, note that if we describe the boundary curve as a + // polynomial in tangential coordinate @p{t=0..1} (along the line) + // and @p{s} in normal direction, then the cubic mapping is such + // that @p{s = a*t**3 + b*t**2 + c*t + d}, and we want to determine + // the interpolation points at @p{t=0.276} and @p{t=0.724} + // (Gauss-Lobatto points). Since at @p{t=0,1} we want a vertex which + // is actually at the boundary, we know that @p{d=0} and @p{a=-b-c}, + // which gives @p{s(0.276)} and @{s(0.726)} in terms of @p{b,c}. As + // side-conditions, we want that the derivatives at @p{t=0} and + // @p{t=1}, i.e. at the vertices match those returned by the + // boundary. // - // The task is then first - // to determine the - // coefficients from the - // tangentials. for that, - // first rotate the - // tangents of @p{s(t)} - // into the global - // coordinate system. they - // are @p{A (1,c)} and @p{A - // (1,-b-2c)} with @p{A} the - // rotation matrix, since - // the tangentials in the - // coordinate system - // relative to the line are - // @p{(1,c)} and @p{(1,-b-2c)} - // at the two vertices, - // respectively. We then - // have to make sure by - // matching @p{b,c} that - // these tangentials are - // orthogonal to the normals - // returned by the boundary - // object + // The task is then first to determine the coefficients from the + // tangentials. for that, first rotate the tangents of @p{s(t)} into + // the global coordinate system. they are @p{A (1,c)} and @p{A + // (1,-b-2c)} with @p{A} the rotation matrix, since the tangentials + // in the coordinate system relative to the line are @p{(1,c)} and + // @p{(1,-b-2c)} at the two vertices, respectively. We then have to + // make sure by matching @p{b,c} that these tangentials are + // orthogonal to the normals returned by the boundary object const Tensor<1,2> coordinate_vector = line->vertex(1) - line->vertex(0); const double h = std::sqrt(coordinate_vector * coordinate_vector); Tensor<1,2> coordinate_axis = coordinate_vector; @@ -151,18 +109,18 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce -face_vertex_normals[1][0] * std::sin(alpha))) -2*c; + QGaussLobatto<1> quad_points(4); + const double t1 = quad_points.point(1)[0]; + const double t2 = quad_points.point(2)[0]; + const double s_t1 = (((-b-c)*t1+b)*t1+c)*t1; + const double s_t2 = (((-b-c)*t2+b)*t2+c)*t2; - // next evaluate the so - // determined cubic - // polynomial at the points - // 1/3 and 2/3, first in - // unit coordinates - const Point<2> new_unit_points[2] = { Point<2>(1./3., 1./27.*(2*b+8*c)), - Point<2>(2./3., 4./27.*b+10./27.*c) + // next evaluate the so determined cubic polynomial at the points + // 1/3 and 2/3, first in unit coordinates + const Point<2> new_unit_points[2] = { Point<2>(t1, s_t1), + Point<2>(t2, s_t2) }; - // then transform these - // points to real - // coordinates by rotating, + // then transform these points to real coordinates by rotating, // scaling and shifting for (unsigned int i=0; i<2; ++i) { diff --git a/deal.II/source/fe/mapping_q.cc b/deal.II/source/fe/mapping_q.cc index 3269194d14..4ee1d4374c 100644 --- a/deal.II/source/fe/mapping_q.cc +++ b/deal.II/source/fe/mapping_q.cc @@ -134,7 +134,7 @@ MappingQ::MappingQ (const unsigned int p, // Qp mapping of cells at the boundary. const QGaussLobatto<1> points(degree+1); tensor_pols = new TensorProductPolynomials - (Polynomials::LagrangeEquidistant::generate_complete_basis(degree)); + (Polynomials::generate_complete_Lagrange_basis(points.get_points())); Assert (n_shape_functions==tensor_pols->n(), ExcInternalError()); Assert(n_inner+n_outer==n_shape_functions, ExcInternalError()); @@ -290,10 +290,8 @@ MappingQ::get_subface_data (const UpdateFlags update_flags, } -// Note that the CellSimilarity flag is -// modifyable, since MappingQ can need to -// recalculate data even when cells are -// similar. +// Note that the CellSimilarity flag is modifyable, since MappingQ can need to +// recalculate data even when cells are similar. template void MappingQ::fill_fe_values ( @@ -308,35 +306,23 @@ MappingQ::fill_fe_values ( std::vector > &normal_vectors, CellSimilarity::Similarity &cell_similarity) const { - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible + // convert data object to internal data for this class. fails with an + // exception if that is not possible Assert (dynamic_cast (&mapping_data) != 0, ExcInternalError()); InternalData &data = static_cast (mapping_data); - // check whether this cell needs - // the full mapping or can be - // treated by a reduced Q1 mapping, - // e.g. if the cell is in the - // interior of the domain + // check whether this cell needs the full mapping or can be treated by a + // reduced Q1 mapping, e.g. if the cell is in the interior of the domain data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - // depending on this result, use this or - // the other data object for the - // mapping. furthermore, we need to - // ensure that the flag indicating - // whether we can use some similarity has - // to be modified - for a general - // MappingQ, the data needs to be - // recomputed anyway since then the - // mapping changes the data. this needs - // to be known also for later operations, - // so modify the variable here. this also - // affects the calculation of the next - // cell -- if we use Q1 data on the next - // cell, the data will still be invalid. + // depending on this result, use this or the other data object for the + // mapping. furthermore, we need to ensure that the flag indicating whether + // we can use some similarity has to be modified - for a general MappingQ, + // the data needs to be recomputed anyway since then the mapping changes the + // data. this needs to be known also for later operations, so modify the + // variable here. this also affects the calculation of the next cell -- if + // we use Q1 data on the next cell, the data will still be invalid. typename MappingQ1::InternalData *p_data=0; if (data.use_mapping_q1_on_current_cell) p_data=&data.mapping_q1_data; @@ -368,34 +354,23 @@ MappingQ::fill_fe_face_values ( std::vector > &exterior_forms, std::vector > &normal_vectors) const { - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible + // convert data object to internal data for this class. fails with an + // exception if that is not possible Assert (dynamic_cast (&mapping_data) != 0, ExcInternalError()); InternalData &data = static_cast (mapping_data); - // check whether this cell needs - // the full mapping or can be - // treated by a reduced Q1 mapping, - // e.g. if the cell is entirely in - // the interior of the domain. note - // that it is not sufficient to ask - // whether the present _face_ is in - // the interior, as the mapping on - // the face depends on the mapping - // of the cell, which in turn - // depends on the fact whether - // _any_ of the faces of this cell - // is at the boundary, not only the - // present face + // check whether this cell needs the full mapping or can be treated by a + // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the + // domain. note that it is not sufficient to ask whether the present _face_ + // is in the interior, as the mapping on the face depends on the mapping of + // the cell, which in turn depends on the fact whether _any_ of the faces of + // this cell is at the boundary, not only the present face data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - // depending on this result, use - // this or the other data object - // for the mapping + // depending on this result, use this or the other data object for the + // mapping typename MappingQ1::InternalData *p_data=0; if (data.use_mapping_q1_on_current_cell) p_data=&data.mapping_q1_data; @@ -430,34 +405,23 @@ MappingQ::fill_fe_subface_values (const typename Triangulation > &exterior_forms, std::vector > &normal_vectors) const { - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible + // convert data object to internal data for this class. fails with an + // exception if that is not possible Assert (dynamic_cast (&mapping_data) != 0, ExcInternalError()); InternalData &data = static_cast (mapping_data); - // check whether this cell needs - // the full mapping or can be - // treated by a reduced Q1 mapping, - // e.g. if the cell is entirely in - // the interior of the domain. note - // that it is not sufficient to ask - // whether the present _face_ is in - // the interior, as the mapping on - // the face depends on the mapping - // of the cell, which in turn - // depends on the fact whether - // _any_ of the faces of this cell - // is at the boundary, not only the - // present face + // check whether this cell needs the full mapping or can be treated by a + // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the + // domain. note that it is not sufficient to ask whether the present _face_ + // is in the interior, as the mapping on the face depends on the mapping of + // the cell, which in turn depends on the fact whether _any_ of the faces of + // this cell is at the boundary, not only the present face data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - // depending on this result, use - // this or the other data object - // for the mapping + // depending on this result, use this or the other data object for the + // mapping typename MappingQ1::InternalData *p_data=0; if (data.use_mapping_q1_on_current_cell) p_data=&data.mapping_q1_data; @@ -499,25 +463,19 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const const unsigned int n_inner_2d=(degree-1)*(degree-1); const unsigned int n_outer_2d=4+4*(degree-1); - // first check whether we have precomputed - // the values for some polynomial degree; - // the sizes of arrays is - // n_inner_2d*n_outer_2d + // first check whether we have precomputed the values for some polynomial + // degree; the sizes of arrays is n_inner_2d*n_outer_2d double const *loqv_ptr=0; switch (degree) { - // for degree==1, we shouldn't have to - // compute any support points, since + // for degree==1, we shouldn't have to compute any support points, since // all of them are on the vertices case 2: { - // (checked these values against the - // output of compute_laplace_vector - // again, and found they're indeed - // right -- just in case someone - // wonders where they come from -- - // WB) + // (checked these values against the output of compute_laplace_vector + // again, and found they're indeed right -- just in case someone wonders + // where they come from -- WB) static const double loqv2[1*8] = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.}; loqv_ptr=&loqv2[0]; @@ -528,133 +486,14 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const break; } - case 3: - { - // (same as above) - static const double loqv3[4*12] - = {80/1053., 1/81., 1/81., 11/1053., 25/117., 44/351., - 7/117., 16/351., 25/117., 44/351., 7/117., 16/351., - 1/81., 80/1053., 11/1053., 1/81., 7/117., 16/351., - 25/117., 44/351., 44/351., 25/117., 16/351., 7/117., - 1/81., 11/1053., 80/1053., 1/81., 44/351., 25/117., - 16/351., 7/117., 7/117., 16/351., 25/117., 44/351., - 11/1053., 1/81., 1/81., 80/1053., 16/351., 7/117., - 44/351., 25/117., 16/351., 7/117., 44/351., 25/117. - }; - Assert (sizeof(loqv3)/sizeof(loqv3[0]) == - n_inner_2d * n_outer_2d, - ExcInternalError()); - - loqv_ptr=&loqv3[0]; - - break; - } - - - case 4: - { - static const double loqv4[9*16] - = {0.07405921850311571, -0.001075744628905992, - -0.001075744628906007, 0.001914292239071463, - 0.2231273865431892, 0.1346851306015187, - 0.03812914216116724, 0.02913160002633252, - 0.02200737428129396, 0.01600835564431224, - 0.2231273865431891, 0.1346851306015187, - 0.03812914216116723, 0.02913160002633253, - 0.02200737428129391, 0.01600835564431222, - - 0.00664803151334206, 0.006648031513342719, - 0.002873452861657458, 0.002873452861657626, - 0.07903572682584378, 0.05969238281250031, - 0.03619864817415824, 0.07903572682584187, - 0.0596923828124999, 0.03619864817415815, - 0.1527716818820237, 0.2348152760709273, - 0.152771681882024, 0.02496269311797778, - 0.04081948955407129, 0.02496269311797789, - - -0.001075744628906923, 0.07405921850311589, - 0.001914292239071339, -0.001075744628905884, - 0.02913160002633509, 0.02200737428129395, - 0.01600835564431229, 0.2231273865431878, - 0.1346851306015183, 0.0381291421611672, - 0.03812914216116729, 0.1346851306015185, - 0.2231273865431898, 0.01600835564431217, - 0.02200737428129394, 0.02913160002633262, - - 0.006648031513342073, 0.002873452861657473, - 0.006648031513342726, 0.002873452861657636, - 0.1527716818820238, 0.2348152760709273, - 0.152771681882024, 0.02496269311797779, - 0.04081948955407131, 0.0249626931179779, - 0.07903572682584376, 0.05969238281250026, - 0.03619864817415824, 0.07903572682584187, - 0.0596923828124998, 0.0361986481741581, - - 0.01106770833333302, 0.01106770833333336, - 0.01106770833333337, 0.01106770833333374, - 0.06770833333333424, 0.1035156250000011, - 0.0677083333333344, 0.06770833333333376, - 0.103515624999999, 0.06770833333333399, - 0.06770833333333422, 0.1035156250000009, - 0.06770833333333436, 0.0677083333333337, - 0.1035156249999988, 0.0677083333333339, - - 0.002873452861657185, 0.006648031513342362, - 0.002873452861657334, 0.006648031513343038, - 0.02496269311797779, 0.04081948955407401, - 0.02496269311797788, 0.1527716818820234, - 0.234815276070926, 0.1527716818820237, - 0.03619864817415819, 0.05969238281250028, - 0.07903572682584407, 0.03619864817415804, - 0.05969238281249986, 0.0790357268258422, - - -0.001075744628906913, 0.00191429223907134, - 0.07405921850311592, -0.001075744628905865, - 0.03812914216116729, 0.1346851306015185, - 0.2231273865431899, 0.01600835564431217, - 0.02200737428129396, 0.02913160002633264, - 0.02913160002633509, 0.02200737428129391, - 0.01600835564431228, 0.2231273865431878, - 0.1346851306015183, 0.03812914216116718, - - 0.002873452861657176, 0.002873452861657321, - 0.006648031513342374, 0.006648031513343037, - 0.03619864817415817, 0.05969238281250032, - 0.07903572682584409, 0.03619864817415805, - 0.05969238281249992, 0.07903572682584221, - 0.02496269311797776, 0.04081948955407392, - 0.02496269311797785, 0.1527716818820233, - 0.2348152760709258, 0.1527716818820236, - - 0.001914292239071237, -0.001075744628906803, - -0.001075744628906778, 0.07405921850311617, - 0.01600835564431228, 0.02200737428129401, - 0.02913160002633524, 0.03812914216116726, - 0.1346851306015182, 0.2231273865431886, - 0.01600835564431228, 0.02200737428129397, - 0.02913160002633523, 0.03812914216116726, - 0.1346851306015181, 0.2231273865431886, - }; - - Assert (sizeof(loqv4)/sizeof(loqv4[0]) == - n_inner_2d * n_outer_2d, - ExcInternalError()); - - loqv_ptr=&loqv4[0]; - - break; - } - - // no other cases implemented, - // so simply fall through + // no other cases implemented, so simply fall through default: break; } if (loqv_ptr!=0) { - // precomputed. copy values to - // the loqvs array + // precomputed. copy values to the loqvs array loqvs.reinit(n_inner_2d, n_outer_2d); for (unsigned int unit_point=0; unit_point::set_laplace_on_quad_vector(Table<2,double> &loqvs) const else { // not precomputed, then do so now - if (dim==2) + if (dim == 2) compute_laplace_vector(loqvs); - else - // computing the Laplace vector for - // faces is not supported in 3d at - // present. presumably, doing so - // would not be so hard: we would - // only have to call the function in - // 2d, i.e. the quad(=face) values in - // 3d are equal to the quad(=cell) - // values in 2d. however, that would - // require us to link in the 2d - // library, which is kind of awkward - // (note that compute_laplace_vector - // really makes use of a lot of 2d - // stuff, such as FEValues etc). an - // alternative would be to precompute - // the values of this array for a - // couple of higher mapping orders, - // pin down their values and insert - // them into the array above. - Assert (false, ExcNotImplemented()); + else if (dim == 3) + { + MappingQ<2,2> mapping_2d(this->degree); + loqvs = mapping_2d.laplace_on_quad_vector; + } } - // the sum of weights of the points - // at the outer rim should be - // one. check this + // the sum of weights of the points at the outer rim should be one. check + // this for (unsigned int unit_point=0; unit_pointdegree, @@ -704,9 +527,8 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const { Assert(degree>1, ExcInternalError()); - // first check whether we have - // precomputed the values for some - // polynomial degree + // first check whether we have precomputed the values for some polynomial + // degree double const *lohv_ptr=0; if (degree==2) { @@ -722,8 +544,7 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const if (lohv_ptr!=0) { - // precomputed. copy values to - // the lohvs array + // precomputed. copy values to the lohvs array lohvs.reinit(n_inner, n_outer); for (unsigned int unit_point=0; unit_point::set_laplace_on_hex_vector(Table<2,double> &lohvs) const // not precomputed, then do so now compute_laplace_vector(lohvs); - // the sum of weights of the points - // at the outer rim should be - // one. check this + // the sum of weights of the points at the outer rim should be one. check + // this for (unsigned int unit_point=0; unit_pointdegree*this->degree, @@ -769,14 +589,11 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const Assert(lvs.n_rows()==0, ExcInternalError()); Assert(dim==2 || dim==3, ExcNotImplemented()); - // for degree==1, we shouldn't have to - // compute any support points, since all of - // them are on the vertices + // for degree==1, we shouldn't have to compute any support points, since all + // of them are on the vertices Assert(degree>1, ExcInternalError()); - // compute the shape - // gradients at the quadrature - // points on the unit cell + // compute the shape gradients at the quadrature points on the unit cell const QGauss quadrature(degree+1); const unsigned int n_q_points=quadrature.size(); @@ -784,8 +601,7 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points); this->compute_shapes(quadrature.get_points(), quadrature_data); - // Compute the stiffness matrix of - // the inner dofs + // Compute the stiffness matrix of the inner dofs FullMatrix S(n_inner); for (unsigned int point=0; point::compute_laplace_vector(Table<2,double> &lvs) const S(i,j) += res * (long double)quadrature.weight(point); } - // Compute the components of T to be the - // product of gradients of inner and + // Compute the components of T to be the product of gradients of inner and // outer shape functions. FullMatrix T(n_inner, n_outer); for (unsigned int point=0; point::compute_laplace_vector(Table<2,double> &lvs) const // S:=S_1*T S_1.mmult(S_1_T,T); - // Resize and initialize the - // lvs + // Resize and initialize the lvs lvs.reinit (n_inner, n_outer); for (unsigned int i=0; i::apply_laplace_vector(const Table<2,double> &lvs, std::vector > &a) const { - // check whether the data we need - // is really available. if you fail - // here and if - // lvs==laplace_on_quad_vector in - // the calling function, then we - // didn't compute the quad laplace - // vector. this is mentioned in the - // constructor of this class, - // although I don't understand the - // reason for not aborting there + // check whether the data we need is really available. if you fail here and + // if lvs==laplace_on_quad_vector in the calling function, then we didn't + // compute the quad laplace vector. this is mentioned in the constructor of + // this class, although I don't understand the reason for not aborting there // any more [WB] Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(degree)); @@ -858,12 +666,9 @@ MappingQ::apply_laplace_vector(const Table<2,double> &lvs, Assert(a.size()==n_outer_apply, ExcDimensionMismatch(a.size(), n_outer_apply)); - // compute each inner point as - // linear combination of the outer - // points. the weights are given by - // the lvs entries, the outer - // points are the first (existing) - // elements of a + // compute each inner point as linear combination of the outer points. the + // weights are given by the lvs entries, the outer points are the first + // (existing) elements of a for (unsigned int unit_point=0; unit_point::compute_mapping_support_points( const typename Triangulation::cell_iterator &cell, std::vector > &a) const { - // if this is a cell for which we - // want to compute the full - // mapping, then get them from the - // following function + // if this is a cell for which we want to compute the full mapping, then get + // them from the following function if (use_mapping_q_on_all_cells || cell->has_boundary_lines()) compute_support_points_laplace(cell, a); else - // otherwise: use a Q1 mapping - // for which the mapping shape - // function support points are - // simply the vertices of the - // cell + // otherwise: use a Q1 mapping for which the mapping shape function + // support points are simply the vertices of the cell { a.resize(GeometryInfo::vertices_per_cell); @@ -908,8 +708,7 @@ void MappingQ::compute_support_points_laplace(const typename Triangulation::cell_iterator &cell, std::vector > &a) const { - // in any case, we need the - // vertices first + // in any case, we need the vertices first a.resize(GeometryInfo::vertices_per_cell); for (unsigned int i=0; i::vertices_per_cell; ++i) a[i] = cell->vertex(i); @@ -921,11 +720,8 @@ MappingQ::compute_support_points_laplace(const typename Triangulat add_line_support_points(cell, a); break; case 2: - // in 2d, add the - // points on the four - // bounding lines to - // the exterior (outer) - // points + // in 2d, add the points on the four bounding lines to the exterior + // (outer) points add_line_support_points (cell, a); if (dim != spacedim) add_quad_support_points(cell, a); @@ -935,9 +731,7 @@ MappingQ::compute_support_points_laplace(const typename Triangulat case 3: { - // in 3d also add the - // points located on - // the boundary faces + // in 3d also add the points located on the boundary faces add_line_support_points (cell, a); add_quad_support_points (cell, a); apply_laplace_vector (laplace_on_hex_vector, a); @@ -956,8 +750,7 @@ void MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &, std::vector > &) const { - // there are no points on bounding - // lines which are to be added + // there are no points on bounding lines which are to be added const unsigned int dim=1; Assert (dim > 1, ExcImpossibleInDim(dim)); } @@ -971,8 +764,7 @@ MappingQ<1,2>::add_line_support_points (const Triangulation<1,2>::cell_iterator { const unsigned int dim = 1; const unsigned int spacedim = 2; - // Ask for the mid point, if that's - // the only thing we need. + // Ask for the mid point, if that's the only thing we need. if (degree==2) { const Boundary *const boundary @@ -980,10 +772,8 @@ MappingQ<1,2>::add_line_support_points (const Triangulation<1,2>::cell_iterator a.push_back(boundary->get_new_point_on_line(cell)); } else - // otherwise call the more - // complicated functions and ask - // for inner points from the - // boundary description + // otherwise call the more complicated functions and ask for inner points + // from the boundary description { std::vector > line_points (degree-1); @@ -1005,8 +795,7 @@ MappingQ<1,3>::add_line_support_points (const Triangulation<1,3>::cell_iterator { const unsigned int dim = 1; const unsigned int spacedim = 3; - // Ask for the mid point, if that's - // the only thing we need. + // Ask for the mid point, if that's the only thing we need. if (degree==2) { const Boundary *const boundary @@ -1014,10 +803,8 @@ MappingQ<1,3>::add_line_support_points (const Triangulation<1,3>::cell_iterator a.push_back(boundary->get_new_point_on_line(cell)); } else - // otherwise call the more - // complicated functions and ask - // for inner points from the - // boundary description + // otherwise call the more complicated functions and ask for inner points + // from the boundary description { std::vector > line_points (degree-1); @@ -1038,8 +825,7 @@ MappingQ::add_line_support_points (const typename Triangulation > &a) const { static const StraightBoundary straight_boundary; - // if we only need the midpoint, - // then ask for it. + // if we only need the midpoint, then ask for it. if (degree==2) { for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) @@ -1056,19 +842,13 @@ MappingQ::add_line_support_points (const typename Triangulation > line_points (degree-1); - // loop over each of the lines, - // and if it is at the - // boundary, then first get the - // boundary description and - // second compute the points on - // it + // loop over each of the lines, and if it is at the boundary, then first + // get the boundary description and second compute the points on it for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) { const typename Triangulation::line_iterator line = cell->line(line_no); @@ -1083,8 +863,7 @@ MappingQ::add_line_support_points (const typename Triangulationget_intermediate_points_on_line (line, line_points); if (dim==3) { - // in 3D, lines might be in wrong - // orientation. if so, reverse + // in 3D, lines might be in wrong orientation. if so, reverse // the vector if (cell->line_orientation(line_no)) a.insert (a.end(), line_points.begin(), line_points.end()); @@ -1092,9 +871,8 @@ MappingQ::add_line_support_points (const typename Triangulation::cell_iterator &cell, vertices_per_cell = GeometryInfo<3>::vertices_per_cell; static const StraightBoundary<3> straight_boundary; - // used if face quad at boundary or - // entirely in the interior of the - // domain + // used if face quad at boundary or entirely in the interior of the domain std::vector > quad_points ((degree-1)*(degree-1)); - // used if only one line of face - // quad is at boundary + // used if only one line of face quad is at boundary std::vector > b(4*degree); - // loop over all faces and collect - // points on them + // loop over all faces and collect points on them for (unsigned int face_no=0; face_no::face_iterator face = cell->face(face_no); - // select the correct mappings - // for the present face + // select the correct mappings for the present face const bool face_orientation = cell->face_orientation(face_no), face_flip = cell->face_flip (face_no), face_rotation = cell->face_rotation (face_no); @@ -1146,9 +919,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, face_rotation)), ExcInternalError()); - // indices of the lines that - // bound a face are given by - // GeometryInfo<3>:: + // indices of the lines that bound a face are given by GeometryInfo<3>:: // face_to_cell_lines for (unsigned int i=0; iline(i)==cell->line(GeometryInfo<3>::face_to_cell_lines( @@ -1156,22 +927,16 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, ExcInternalError()); #endif - // if face at boundary, then - // ask boundary object to - // return intermediate points - // on it + // if face at boundary, then ask boundary object to return intermediate + // points on it if (face->at_boundary()) { face->get_triangulation().get_boundary(face->boundary_indicator()) .get_intermediate_points_on_quad (face, quad_points); - // in 3D, the orientation, flip and - // rotation of the face might not - // match what we expect here, namely - // the standard orientation. thus - // reorder points accordingly. since - // a Mapping uses the same shape - // function as an FEQ, we can ask a - // FEQ to do the reordering for us. + // in 3D, the orientation, flip and rotation of the face might not + // match what we expect here, namely the standard orientation. thus + // reorder points accordingly. since a Mapping uses the same shape + // function as an FEQ, we can ask a FEQ to do the reordering for us. for (unsigned int i=0; i::cell_iterator &cell, } else { - // face is not at boundary, - // but maybe some of its - // lines are. count them + // face is not at boundary, but maybe some of its lines are. count + // them unsigned int lines_at_boundary=0; for (unsigned int i=0; iline(i)->at_boundary()) @@ -1190,37 +954,23 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, Assert(lines_at_boundary<=lines_per_face, ExcInternalError()); - // if at least one of the - // lines bounding this quad - // is at the boundary, then - // collect points - // separately + // if at least one of the lines bounding this quad is at the + // boundary, then collect points separately if (lines_at_boundary>0) { - // call of function - // apply_laplace_vector - // increases size of b - // about 1. There - // resize b for the - // case the mentioned - // function was already - // called. + // call of function apply_laplace_vector increases size of b + // about 1. There resize b for the case the mentioned function + // was already called. b.resize(4*degree); - // b is of size - // 4*degree, make sure - // that this is the - // right size + // b is of size 4*degree, make sure that this is the right size Assert(b.size()==vertices_per_face+lines_per_face*(degree-1), ExcDimensionMismatch(b.size(), vertices_per_face+lines_per_face*(degree-1))); - // sort the points into b. We - // used access from the cell (not - // from the face) to fill b, so - // we can assume a standard face - // orientation. Doing so, the - // calculated points will be in + // sort the points into b. We used access from the cell (not + // from the face) to fill b, so we can assume a standard face + // orientation. Doing so, the calculated points will be in // standard orientation as well. for (unsigned int i=0; i::face_to_cell_vertices(face_no, i)]; @@ -1231,8 +981,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, a[vertices_per_cell + GeometryInfo<3>::face_to_cell_lines( face_no, i)*(degree-1)+j]; - // Now b includes the support - // points on the quad and we can + // Now b includes the support points on the quad and we can // apply the laplace vector apply_laplace_vector(laplace_on_quad_vector, b); Assert(b.size()==4*degree+(degree-1)*(degree-1), @@ -1243,22 +992,14 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, } else { - // face is entirely in - // the interior. get - // intermediate points - // from a straight - // boundary object + // face is entirely in the interior. get intermediate points + // from a straight boundary object straight_boundary.get_intermediate_points_on_quad (face, quad_points); - // in 3D, the orientation, flip - // and rotation of the face might - // not match what we expect here, - // namely the standard - // orientation. thus reorder - // points accordingly. since a - // Mapping uses the same shape - // function as an FEQ, we can ask - // a FEQ to do the reordering for - // us. + // in 3D, the orientation, flip and rotation of the face might + // not match what we expect here, namely the standard + // orientation. thus reorder points accordingly. since a Mapping + // uses the same shape function as an FEQ, we can ask a FEQ to + // do the reordering for us. for (unsigned int i=0; i::transform ( const MappingType mapping_type) const { AssertDimension (input.size(), output.size()); - // The data object may be just a - // MappingQ1::InternalData, so we - // have to test for this first. + // The data object may be just a MappingQ1::InternalData, so we have to test + // for this first. const typename MappingQ1::InternalData *q1_data = dynamic_cast::InternalData *> (&mapping_data); Assert(q1_data!=0, ExcInternalError()); - // If it is a genuine - // MappingQ::InternalData, we have - // to test further + // If it is a genuine MappingQ::InternalData, we have to test further if (!q1_data->is_mapping_q1_data) { Assert (dynamic_cast(&mapping_data) != 0, ExcInternalError()); const InternalData &data = static_cast(mapping_data); - // If we only use the - // Q1-portion, we have to - // extract that data object + // If we only use the Q1-portion, we have to extract that data object if (data.use_mapping_q1_on_current_cell) q1_data = &data.mapping_q1_data; } - // Now, q1_data should have the - // right tensors in it and we call - // the base classes transform - // function + // Now, q1_data should have the right tensors in it and we call the base + // classes transform function MappingQ1::transform(input, output, *q1_data, mapping_type); } @@ -1346,31 +1080,24 @@ MappingQ::transform ( const MappingType mapping_type) const { AssertDimension (input.size(), output.size()); - // The data object may be just a - // MappingQ1::InternalData, so we - // have to test for this first. + // The data object may be just a MappingQ1::InternalData, so we have to test + // for this first. const typename MappingQ1::InternalData *q1_data = dynamic_cast::InternalData *> (&mapping_data); Assert(q1_data!=0, ExcInternalError()); - // If it is a genuine - // MappingQ::InternalData, we have - // to test further + // If it is a genuine MappingQ::InternalData, we have to test further if (!q1_data->is_mapping_q1_data) { Assert (dynamic_cast(&mapping_data) != 0, ExcInternalError()); const InternalData &data = static_cast(mapping_data); - // If we only use the - // Q1-portion, we have to - // extract that data object + // If we only use the Q1-portion, we have to extract that data object if (data.use_mapping_q1_on_current_cell) q1_data = &data.mapping_q1_data; } - // Now, q1_data should have the - // right tensors in it and we call - // the base classes transform - // function + // Now, q1_data should have the right tensors in it and we call the base + // classes transform function MappingQ1::transform(input, output, *q1_data, mapping_type); } @@ -1383,31 +1110,24 @@ void MappingQ::transform const MappingType mapping_type) const { AssertDimension (input.size(), output.size()); - // The data object may be just a - // MappingQ1::InternalData, so we - // have to test for this first. + // The data object may be just a MappingQ1::InternalData, so we have to test + // for this first. const typename MappingQ1::InternalData *q1_data = dynamic_cast::InternalData *> (&mapping_data); Assert(q1_data!=0, ExcInternalError()); - // If it is a genuine - // MappingQ::InternalData, we have - // to test further + // If it is a genuine MappingQ::InternalData, we have to test further if (!q1_data->is_mapping_q1_data) { Assert (dynamic_cast(&mapping_data) != 0, ExcInternalError()); const InternalData &data = static_cast(mapping_data); - // If we only use the - // Q1-portion, we have to - // extract that data object + // If we only use the Q1-portion, we have to extract that data object if (data.use_mapping_q1_on_current_cell) q1_data = &data.mapping_q1_data; } - // Now, q1_data should have the - // right tensors in it and we call - // the base classes transform - // function + // Now, q1_data should have the right tensors in it and we call the base + // classes transform function MappingQ1::transform(input, output, *q1_data, mapping_type); } @@ -1418,11 +1138,9 @@ MappingQ:: transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, const Point &p) const { - // Use the get_data function to - // create an InternalData with data - // vectors of the right size and - // transformation shape values - // already computed at point p. + // Use the get_data function to create an InternalData with data vectors of + // the right size and transformation shape values already computed at point + // p. const Quadrature point_quadrature(p); std::auto_ptr mdata (dynamic_cast ( @@ -1437,8 +1155,7 @@ transform_unit_to_real_cell (const typename Triangulation::cell_it &*mdata); compute_mapping_support_points(cell, p_data->mapping_support_points); - // If this should be Q1, ignore all - // other support points. + // If this should be Q1, ignore all other support points. if (p_data->shape_values.size()mapping_support_points.size()) p_data->mapping_support_points.resize (GeometryInfo::vertices_per_cell); @@ -1455,24 +1172,15 @@ MappingQ:: transform_real_to_unit_cell (const typename Triangulation::cell_iterator &cell, const Point &p) const { - // first a Newton iteration based - // on a Q1 mapping to get a good - // starting point, the idea being - // that this is cheaper than trying - // to start with the real mapping - // and likely also more robust. + // first a Newton iteration based on a Q1 mapping to get a good starting + // point, the idea being that this is cheaper than trying to start with the + // real mapping and likely also more robust. // - // that said, this doesn't always - // work: there are cases where the - // point is outside the cell and - // the inverse mapping doesn't - // converge. in that case, use the - // center point of the cell as a - // starting point if we are to go - // on using the full mapping, or - // just propagate up the exception - // if we had no intention of - // continuing with the full mapping + // that said, this doesn't always work: there are cases where the point is + // outside the cell and the inverse mapping doesn't converge. in that case, + // use the center point of the cell as a starting point if we are to go on + // using the full mapping, or just propagate up the exception if we had no + // intention of continuing with the full mapping Point initial_p_unit; try { @@ -1481,11 +1189,8 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it } catch (const typename Mapping::ExcTransformationFailed &) { - // mirror the conditions of the - // code below to determine if - // we need to use an arbitrary - // starting point or if we just - // need to rethrow the + // mirror the conditions of the code below to determine if we need to + // use an arbitrary starting point or if we just need to rethrow the // exception if (cell->has_boundary_lines() || @@ -1500,30 +1205,22 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it throw; } - // then a Newton iteration based on the - // full MappingQ if we need this. note that - // for interior cells with dim==spacedim, - // the mapping used is in fact a Q1 - // mapping, so there is nothing we need to - // do unless the iteration above failed + // then a Newton iteration based on the full MappingQ if we need this. note + // that for interior cells with dim==spacedim, the mapping used is in fact a + // Q1 mapping, so there is nothing we need to do unless the iteration above + // failed if (cell->has_boundary_lines() || use_mapping_q_on_all_cells || (dim!=spacedim) ) { - // use the full mapping. in case the - // function above should have given us - // something back that lies outside the - // unit cell (that might happen because - // we may have given a point 'p' that - // lies inside the cell with the higher - // order mapping, but outside the - // Q1-mapped reference cell), then - // project it back into the reference - // cell in hopes that this gives a - // better starting point to the - // following iteration + // use the full mapping. in case the function above should have given us + // something back that lies outside the unit cell (that might happen + // because we may have given a point 'p' that lies inside the cell with + // the higher order mapping, but outside the Q1-mapped reference cell), + // then project it back into the reference cell in hopes that this gives + // a better starting point to the following iteration initial_p_unit = GeometryInfo::project_to_unit_cell(initial_p_unit); const Quadrature point_quadrature(initial_p_unit); @@ -1539,13 +1236,11 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it compute_mapping_support_points (cell, mdata->mapping_support_points); - // If this is a q1 mapping, - // then only use the support - // points on the vertices. + // If this is a q1 mapping, then only use the support points on the + // vertices. if (mdata->shape_values.size() < mdata->mapping_support_points.size()) mdata->mapping_support_points.resize(GeometryInfo::vertices_per_cell); - return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit, *mdata); } else diff --git a/deal.II/source/fe/mapping_q_eulerian.cc b/deal.II/source/fe/mapping_q_eulerian.cc index c8ba737474..902c6cd359 100644 --- a/deal.II/source/fe/mapping_q_eulerian.cc +++ b/deal.II/source/fe/mapping_q_eulerian.cc @@ -71,20 +71,15 @@ SupportQuadrature (const unsigned int map_degree) : Quadrature(Utilities::fixed_power(map_degree+1)) { - // first we determine the support points - // on the unit cell in lexicographic order. - // for this purpose we can use an interated - // trapezoidal quadrature rule. - const QTrapez<1> q1d; - const QIterated q_iterated(q1d,map_degree); + // first we determine the support points on the unit cell in lexicographic + // order, which are (in accordance with MappingQ) the support points of + // QGaussLobatto. + const QGaussLobatto q_iterated(map_degree+1); const unsigned int n_q_points = q_iterated.size(); - // we then need to define a renumbering - // vector that allows us to go from a - // lexicographic numbering scheme to a hierarchic - // one. this fragment is taking almost verbatim - // from the MappingQ class. - + // we then need to define a renumbering vector that allows us to go from a + // lexicographic numbering scheme to a hierarchic one. this fragment is + // taking almost verbatim from the MappingQ class. std::vector renumber(n_q_points); std::vector dpo(dim+1, 1U); for (unsigned int i=1; i (dpo, 1, map_degree), renumber); - // finally we assign the quadrature points in the - // required order. - + // finally we assign the quadrature points in the required order. for (unsigned int q=0; qquadrature_points[renumber[q]] = q_iterated.point(q); } @@ -112,43 +105,30 @@ compute_mapping_support_points std::vector > &a) const { - // first, basic assertion - // with respect to vector size, + // first, basic assertion with respect to vector size, const types::global_dof_index n_dofs = euler_dof_handler->n_dofs(); const types::global_dof_index vector_size = euler_vector->size(); AssertDimension(vector_size,n_dofs); - // we then transform our tria iterator - // into a dof iterator so we can - // access data not associated with - // triangulations + // we then transform our tria iterator into a dof iterator so we can access + // data not associated with triangulations typename DoFHandler::cell_iterator dof_cell(*cell, euler_dof_handler); Assert (dof_cell->active() == true, ExcInactiveCell()); - // our quadrature rule is chosen - // so that each quadrature point - // corresponds to a support point - // in the undeformed configuration. - // we can then query the given - // displacement field at these points - // to determine the shift vector that - // maps the support points to the - // deformed configuration. - - // we assume that the given field contains - // dim displacement components, but - // that there may be other solution - // components as well (e.g. pressures). - // this class therefore assumes that the - // first dim components represent the - // actual shift vector we need, and simply - // ignores any components after that. - // this implies that the user should order - // components appropriately, or create a - // separate dof handler for the displacements. + // our quadrature rule is chosen so that each quadrature point corresponds + // to a support point in the undeformed configuration. we can then query + // the given displacement field at these points to determine the shift + // vector that maps the support points to the deformed configuration. + + // we assume that the given field contains dim displacement components, but + // that there may be other solution components as well (e.g. pressures). + // this class therefore assumes that the first dim components represent the + // actual shift vector we need, and simply ignores any components after + // that. this implies that the user should order components appropriately, + // or create a separate dof handler for the displacements. const unsigned int n_support_pts = support_quadrature.size(); const unsigned int n_components = euler_dof_handler->get_fe().n_components(); @@ -157,20 +137,15 @@ compute_mapping_support_points std::vector > shift_vector(n_support_pts,Vector(n_components)); - // fill shift vector for each - // support point using an fe_values - // object. make sure that the - // fe_values variable isn't used - // simulatenously from different + // fill shift vector for each support point using an fe_values object. make + // sure that the fe_values variable isn't used simulatenously from different // threads Threads::Mutex::ScopedLock lock(fe_values_mutex); fe_values.reinit(dof_cell); fe_values.get_function_values(*euler_vector, shift_vector); - // and finally compute the positions of the - // support points in the deformed + // and finally compute the positions of the support points in the deformed // configuration. - a.resize(n_support_pts); for (unsigned int q=0; q::fill_fe_values ( std::vector > &normal_vectors, CellSimilarity::Similarity &cell_similarity) const { - // disable any previously detected - // similarity and hand on to the respective + // disable any previously detected similarity and hand on to the respective // function of the base class. cell_similarity = CellSimilarity::invalid_next_cell; MappingQ::fill_fe_values (cell, q, mapping_data, diff --git a/deal.II/source/grid/tria_boundary.cc b/deal.II/source/grid/tria_boundary.cc index f956c1d16c..65b3bafe7b 100644 --- a/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/source/grid/tria_boundary.cc @@ -247,6 +247,33 @@ project_to_surface (const typename Triangulation::hex_iterator &, } + +template +const std::vector > & +Boundary:: +get_line_support_points (const unsigned int n_intermediate_points) const +{ + if (points.size() <= n_intermediate_points || + points[n_intermediate_points].get() == 0) + { + Threads::Mutex::ScopedLock lock(mutex); + if (points.size() <= n_intermediate_points) + points.resize(n_intermediate_points+1); + + // another thread might have created points in the meantime + if (points[n_intermediate_points].get() == 0) + { + std_cxx1x::shared_ptr > + quadrature (new QGaussLobatto<1>(n_intermediate_points+2)); + points[n_intermediate_points] = quadrature; + } + } + return points[n_intermediate_points]->get_points(); +} + + + + /* -------------------------- StraightBoundary --------------------- */ @@ -394,15 +421,18 @@ get_intermediate_points_on_line (const Triangulation<1, 2>::line_iterator &line, const unsigned int n=points.size(); Assert(n>0, ExcInternalError()); - const double dx=1./(n+1); - double x=dx; - + // Use interior points of QGaussLobatto quadrature formula support points + // for consistency with MappingQ + const std::vector > &line_points = this->get_line_support_points(n); const Point vertices[2] = { line->vertex(0), line->vertex(1) }; - for (unsigned int i=0; i::li const unsigned int n=points.size(); Assert(n>0, ExcInternalError()); - const double dx=1./(n+1); - double x=dx; + // Use interior points of QGaussLobatto quadrature formula support points + // for consistency with MappingQ + const std::vector > &line_points = this->get_line_support_points(n); const Point vertices[2] = { line->vertex(0), line->vertex(1) }; - for (unsigned int i=0; i::quad_iterator &quad, // is n a square number Assert(m*m==n, ExcInternalError()); - const double ds=1./(m+1); - double y=ds; + const std::vector > &line_points = this->get_line_support_points(m); const Point vertices[4] = { quad->vertex(0), quad->vertex(1), @@ -464,14 +497,17 @@ get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, quad->vertex(3) }; - for (unsigned int i=0; i::quad_iterator &quad, // is n a square number Assert(m*m==n, ExcInternalError()); - const double ds=1./(m+1); - double y=ds; + const std::vector > &line_points = this->get_line_support_points(m); const Point vertices[4] = { quad->vertex(0), quad->vertex(1), @@ -499,14 +534,17 @@ get_intermediate_points_on_quad (const Triangulation<2,3>::quad_iterator &quad, quad->vertex(3) }; - for (unsigned int i=0; i::get_intermediate_points_between_points ( const unsigned int n=points.size(); Assert(n>0, ExcInternalError()); - // Do a simple linear interpolation - // followed by projection, using - // the same algorithm as above - const Point ds = (v1-v0) / (n+1); + // Do a simple linear interpolation followed by projection, using the same + // algorithm as above + const std::vector > &line_points = this->get_line_support_points(n); for (unsigned int i=0; i middle = v0 + (i+1)*ds; + const double x = line_points[i+1][0]; + const Point middle = (1-x)*v0 + x*v1; const Point vector_from_axis = (middle-point_on_axis) - ((middle-point_on_axis) * direction) * direction; @@ -307,24 +307,23 @@ get_intermediate_points_between_points (const Point &p0, { const unsigned int n = points.size (); const Point axis = x_1 - x_0; - const Point dx = (p1 - p0) / (n + 1); Assert (n > 0, ExcInternalError ()); - for (unsigned int i = 0; i < n; ++i) + const std::vector > &line_points = this->get_line_support_points(n); + + for (unsigned int i=0; i x_i = p0 + (i + 1) * dx; - // To project this point on the - // boundary of the cone we first - // compute the orthogonal - // projection of this point onto - // the axis of the cone. + const Point x_i = (1-x)*p0 + x*p1; + // To project this point on the boundary of the cone we first compute + // the orthogonal projection of this point onto the axis of the cone. const double c = (x_i - x_0) * axis / axis.square (); const Point x_ip = x_0 + c * axis; - // Compute the projection of - // the middle point on the - // boundary of the cone. + // Compute the projection of the middle point on the boundary of the + // cone. points[i] = x_ip + get_radius (x_ip) * (x_i - x_ip) / (x_i - x_ip).norm (); } } @@ -337,16 +336,11 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c const Point axis = x_1 - x_0; // Compute the middle point of the line. const Point middle = StraightBoundary::get_new_point_on_line (line); - // To project it on the boundary of - // the cone we first compute the - // orthogonal projection of the - // middle point onto the axis of - // the cone. + // To project it on the boundary of the cone we first compute the orthogonal + // projection of the middle point onto the axis of the cone. const double c = (middle - x_0) * axis / axis.square (); const Point middle_p = x_0 + c * axis; - // Compute the projection of the - // middle point on the boundary - // of the cone. + // Compute the projection of the middle point on the boundary of the cone. return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm (); } @@ -360,20 +354,14 @@ get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const const int dim = 3; const Point axis = x_1 - x_0; - // Compute the middle point of the - // quad. + // Compute the middle point of the quad. const Point middle = StraightBoundary<3>::get_new_point_on_quad (quad); - // Same algorithm as above: To - // project it on the boundary of - // the cone we first compute the - // orthogonal projection of the - // middle point onto the axis of - // the cone. + // Same algorithm as above: To project it on the boundary of the cone we + // first compute the orthogonal projection of the middle point onto the axis + // of the cone. const double c = (middle - x_0) * axis / axis.square (); const Point middle_p = x_0 + c * axis; - // Compute the projection of the - // middle point on the boundary - // of the cone. + // Compute the projection of the middle point on the boundary of the cone. return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm (); } @@ -475,15 +463,12 @@ get_normals_at_vertices (const typename Triangulation::face_iterator &face, for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_cell; ++vertex) { - // Compute the orthogonal - // projection of the vertex onto - // the axis of the cone. + // Compute the orthogonal projection of the vertex onto the axis of the + // cone. const double c = (face->vertex (vertex) - x_0) * axis / axis.square (); const Point vertex_p = x_0 + c * axis; - // Then compute the vector - // pointing from the point - // vertex_p on the axis - // to the vertex. + // Then compute the vector pointing from the point vertex_p on + // the axis to the vertex. const Point axis_to_vertex = face->vertex (vertex) - vertex_p; face_vertex_normals[vertex] = axis_to_vertex / axis_to_vertex.norm (); @@ -629,49 +614,28 @@ HyperBallBoundary::get_intermediate_points_between_points ( Assert(std::fabs(v1.square()-r2) pm=0.5*(v0+v1); const double h=std::sqrt(pm.square()); - double beta=0; - unsigned int left_index=0, right_index=0; - if ((n+1)%2==0) - { - // if the number of - // parts is even insert - // the midpoint - left_index=(n-1)/2; - right_index=left_index; - points[left_index]=pm; - --left_index; - ++right_index; - beta=d_alpha; - } - else - { - beta=(n/2+1)*d_alpha-alpha/2.; - right_index=n/2; - left_index=n/2-1; - } - // n even: m=n/2, // n odd: m=(n-1)/2 + const std::vector > &line_points = this->get_line_support_points(n); const unsigned int m=n/2; - for (unsigned int i=0; i=1, ExcInternalError()); - - double d=h*std::tan(beta); - points[right_index]=pm+d/length*(v1-v0); - points[left_index]=pm-d/length*(v1-v0); + const double beta = alpha * (line_points[i+1][0]-0.5); + const double d=h*std::tan(beta); + points[i]=pm+d/length*(v1-v0); + points[n-1-i]=pm-d/length*(v1-v0); } + if ((n+1)%2==0) + // if the number of parts is even insert the midpoint + points[(n-1)/2]=pm; - // project the points from the - // straight line to the - // HyperBallBoundary + + // project the points from the straight line to the HyperBallBoundary for (unsigned int i=0; i:: get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, std::vector > &points) const { - // check whether center of object is - // at x==0, since then it belongs - // to the plane part of the - // boundary + // check whether center of object is at x==0, since then it belongs to the + // plane part of the boundary const Point line_center = line->center(); if (line_center(0) == this->center(0)) return StraightBoundary::get_intermediate_points_on_line (line, points); @@ -908,10 +870,8 @@ get_intermediate_points_on_quad (const typename Triangulation::quad_iterato points[0]=get_new_point_on_quad(quad); else { - // check whether center of - // object is at x==0, since - // then it belongs to the plane - // part of the boundary + // check whether center of object is at x==0, since then it belongs to + // the plane part of the boundary const Point quad_center = quad->center(); if (quad_center(0) == this->center(0)) StraightBoundary::get_intermediate_points_on_quad (quad, points); @@ -950,10 +910,8 @@ HalfHyperBallBoundary:: get_normals_at_vertices (const typename Triangulation::face_iterator &face, typename Boundary::FaceVertexNormals &face_vertex_normals) const { - // check whether center of object is - // at x==0, since then it belongs - // to the plane part of the - // boundary + // check whether center of object is at x==0, since then it belongs to the + // plane part of the boundary const Point quad_center = face->center(); if (quad_center(0) == this->center(0)) StraightBoundary::get_normals_at_vertices (face, face_vertex_normals); @@ -1005,10 +963,8 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c { switch (dim) { - // in 2d, first check whether the two - // end points of the line are on the - // axis of symmetry. if so, then return - // the mid point + // in 2d, first check whether the two end points of the line are on the + // axis of symmetry. if so, then return the mid point case 2: { if ((line->vertex(0)(0) == this->center(0)) @@ -1016,17 +972,13 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c (line->vertex(1)(0) == this->center(0))) return (line->vertex(0) + line->vertex(1))/2; else - // otherwise we are on the outer or - // inner part of the shell. proceed + // otherwise we are on the outer or inner part of the shell. proceed // as in the base class return HyperShellBoundary::get_new_point_on_line (line); } - // in 3d, a line is a straight - // line if it is on the symmetry - // plane and if not both of its - // end points are on either the - // inner or outer sphere + // in 3d, a line is a straight line if it is on the symmetry plane and if + // not both of its end points are on either the inner or outer sphere case 3: { @@ -1047,8 +999,7 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c - outer_radius) < 1e-12 * outer_radius)))) return (line->vertex(0) + line->vertex(1))/2; else - // otherwise we are on the outer or - // inner part of the shell. proceed + // otherwise we are on the outer or inner part of the shell. proceed // as in the base class return HyperShellBoundary::get_new_point_on_line (line); } @@ -1079,10 +1030,8 @@ Point HalfHyperShellBoundary:: get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const { - // if this quad is on the symmetry plane, - // take the center point and project it - // outward to the same radius as the - // centers of the two radial lines + // if this quad is on the symmetry plane, take the center point and project + // it outward to the same radius as the centers of the two radial lines if ((quad->vertex(0)(0) == this->center(0)) && (quad->vertex(1)(0) == this->center(0)) && (quad->vertex(2)(0) == this->center(0)) && @@ -1119,9 +1068,8 @@ get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) c Assert (false, ExcInternalError()); } - // otherwise we are on the outer or - // inner part of the shell. proceed - // as in the base class + // otherwise we are on the outer or inner part of the shell. proceed as in + // the base class return HyperShellBoundary::get_new_point_on_quad (quad); } @@ -1135,10 +1083,8 @@ get_intermediate_points_on_line (const typename Triangulation::line_iterato { switch (dim) { - // in 2d, first check whether the two - // end points of the line are on the - // axis of symmetry. if so, then return - // the mid point + // in 2d, first check whether the two end points of the line are on the + // axis of symmetry. if so, then return the mid point case 2: { if ((line->vertex(0)(0) == this->center(0)) @@ -1146,18 +1092,14 @@ get_intermediate_points_on_line (const typename Triangulation::line_iterato (line->vertex(1)(0) == this->center(0))) StraightBoundary::get_intermediate_points_on_line (line, points); else - // otherwise we are on the outer or - // inner part of the shell. proceed + // otherwise we are on the outer or inner part of the shell. proceed // as in the base class HyperShellBoundary::get_intermediate_points_on_line (line, points); break; } - // in 3d, a line is a straight - // line if it is on the symmetry - // plane and if not both of its - // end points are on either the - // inner or outer sphere + // in 3d, a line is a straight line if it is on the symmetry plane and if + // not both of its end points are on either the inner or outer sphere case 3: { if (((line->vertex(0)(0) == this->center(0)) @@ -1177,8 +1119,7 @@ get_intermediate_points_on_line (const typename Triangulation::line_iterato - outer_radius) < 1e-12 * outer_radius)))) StraightBoundary::get_intermediate_points_on_line (line, points); else - // otherwise we are on the outer or - // inner part of the shell. proceed + // otherwise we are on the outer or inner part of the shell. proceed // as in the base class HyperShellBoundary::get_intermediate_points_on_line (line, points); @@ -1200,10 +1141,8 @@ get_intermediate_points_on_quad (const typename Triangulation::quad_iterato { Assert (dim < 3, ExcNotImplemented()); - // check whether center of object is - // at x==0, since then it belongs - // to the plane part of the - // boundary + // check whether center of object is at x==0, since then it belongs to the + // plane part of the boundary const Point quad_center = quad->center(); if (quad_center(0) == this->center(0)) StraightBoundary::get_intermediate_points_on_quad (quad, points); @@ -1350,9 +1289,9 @@ TorusBoundary<2,3>::get_new_point_on_line (const Triangulation<2,3>::line_iterat Point<2> middle(0,0); - //Take care for periodic conditions, - //For instance phi0= 0, phi1= 3/2*Pi middle has to be 7/4*Pi not 3/4*Pi - //This also works for -Pi/2 + Pi, middle is 5/4*Pi + //Take care for periodic conditions, For instance phi0= 0, phi1= 3/2*Pi + //middle has to be 7/4*Pi not 3/4*Pi. This also works for -Pi/2 + Pi, middle + //is 5/4*Pi for (unsigned int i=0; i<2; i++) if (std::abs(p0(i)-p1(i))> numbers::PI) middle(i)=2*numbers::PI; @@ -1440,7 +1379,7 @@ TorusBoundary<2,3>:: get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator &line, std::vector< Point< 3 > > &points) const { - //Almost the same implementation as StraightBoundary<2,3> + //Almost the same implementation as StraightBoundary<2,3> unsigned int npoints=points.size(); if (npoints==0) return; @@ -1453,10 +1392,10 @@ get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator &lin offset[0]=0; offset[1]=0; - //Take care for periodic conditions & negative angles, - //see get_new_point_on_line() above - //Because we dont have a symmetric interpolation (just the middle) we need to - //add 2*Pi to each almost zero and negative angles. + //Take care for periodic conditions & negative angles, see + //get_new_point_on_line() above. Because we dont have a symmetric + //interpolation (just the middle) we need to add 2*Pi to each almost zero + //and negative angles. for (unsigned int i=0; i<2; i++) for (unsigned int j=1; j<2; j++) { @@ -1473,11 +1412,11 @@ get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator &lin p[j](i)+=2*numbers::PI*offset[i]; - double dx=1.0/(npoints+1); - double x=dx; Point<2> target; - for (unsigned int i=0; i > &line_points = this->get_line_support_points(npoints); + for (unsigned int i=0; i::quad_iterator &qua // is n a square number Assert(m*m==n, ExcInternalError()); - const double ds=1./(m+1); - double y=ds; - Point<2> p[4]; for (unsigned int i=0; i<4; i++) @@ -1510,10 +1446,10 @@ get_intermediate_points_on_quad (const Triangulation< 2, 3 >::quad_iterator &qua offset[0]=0; offset[1]=0; - //Take care for periodic conditions & negative angles, - //see get_new_point_on_line() above - //Because we dont have a symmetric interpolation (just the middle) we need to - //add 2*Pi to each almost zero and negative angles. + //Take care for periodic conditions & negative angles, see + //get_new_point_on_line() above. Because we dont have a symmetric + //interpolation (just the middle) we need to add 2*Pi to each almost zero + //and negative angles. for (unsigned int i=0; i<2; i++) for (unsigned int j=1; j<4; j++) { @@ -1529,11 +1465,13 @@ get_intermediate_points_on_quad (const Triangulation< 2, 3 >::quad_iterator &qua if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles p[j](i)+=2*numbers::PI*offset[i]; - for (unsigned int i=0; i > &line_points = this->get_line_support_points(m); + for (unsigned int i=0; i