From fd129799847e0372ae021eeb6774cb827a127105 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 13 Mar 2001 12:22:00 +0000 Subject: [PATCH] Minor clean-ups. Prepare Mapping_Q for possible creation of a Mapping_C1 class. git-svn-id: https://svn.dealii.org/trunk@4195 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/mapping.h | 8 +- .../deal.II/include/fe/mapping_cartesian.h | 2 +- deal.II/deal.II/include/fe/mapping_q.h | 79 +++++--- deal.II/deal.II/include/fe/mapping_q1.h | 52 ++++-- deal.II/deal.II/source/fe/mapping_q.cc | 175 ++++++++++++------ deal.II/deal.II/source/fe/mapping_q1.cc | 31 ++-- 6 files changed, 229 insertions(+), 118 deletions(-) diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index dfc0513c51..28c7afbd50 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors +// Copyright (C) 2000, 2001 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -51,6 +51,7 @@ class Mapping : public Subscriptor * Class for internal data of finite * element and mapping objects. */ +//TODO: can we make the following class protected? class InternalDataBase: public Subscriptor { public: @@ -111,7 +112,8 @@ class Mapping : public Subscriptor * Virtual destructor. */ virtual ~Mapping (); - + +//TODO: why make the following functions public? they are only helpful for fevalues and maybe the finite elements? /** * Prepare internal data * structures and fill in values @@ -246,7 +248,7 @@ class Mapping : public Subscriptor virtual void transform_contravariant (std::vector > &dst, const std::vector > &src, const InternalDataBase& internal, - const unsigned int src_offset) const = 0; + const unsigned int src_offset) const = 0; /** * Tranform a field of covariant vectors. diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h index 5b5d82bc13..e32c05fcf4 100644 --- a/deal.II/deal.II/include/fe/mapping_cartesian.h +++ b/deal.II/deal.II/include/fe/mapping_cartesian.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors +// Copyright (C) 2001 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index 4938a56974..88e6ab49ce 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors +// Copyright (C) 2000, 2001 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -57,6 +57,7 @@ class MappingQ : public MappingQ1 */ ~MappingQ (); +//TODO: why make the following functions public? they are only helpful for fevalues and maybe the finite elements? /** * Implementation of the interface in * @ref{Mapping}. @@ -178,7 +179,7 @@ class MappingQ : public MappingQ1 /** * Constructor. */ - InternalData(unsigned int n_shape_functions); + InternalData (const unsigned int n_shape_functions); /** * Unit normal vectors. Used @@ -292,11 +293,18 @@ class MappingQ : public MappingQ1 /** * Takes a - * @p{set_laplace_on_hex(quad)_vector} + * @p{laplace_on_hex(quad)_vector} * and applies it to the vector * @p{a} to compute the inner - * support points. They are - * appended to the vector @p{a}. + * support points as a linear + * combination of the exterior + * points. + * + * The vector @p{a} initially + * containts the locations of the + * @p{n_outer} points, the + * @p{n_inner} computed inner + * points are appended. */ void apply_laplace_vector(const std::vector > &lvs, std::vector > &a) const; @@ -319,9 +327,9 @@ class MappingQ : public MappingQ1 * solution of a Laplace equation * with the position of the outer * support points as boundary - * values. The outer support - * points are all support points - * except of the inner ones. + * values, in order to make the + * transformation as smooth as + * possible. */ void compute_support_points_laplace( const typename Triangulation::cell_iterator &cell, @@ -340,31 +348,59 @@ class MappingQ : public MappingQ1 std::vector > &a) const; /** - * For @p{dim=2,3}. Adds (appends) the - * support points of all lines to - * the vector a. + * For @p{dim=2,3}. Append + * (appends) 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 * @p{compute_support_points_simple(laplace)} * functions. For @p{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. */ - void add_line_support_points (const Triangulation::cell_iterator &cell, - std::vector > &a) const; + virtual void + add_line_support_points (const Triangulation::cell_iterator &cell, + std::vector > &a) const; /** - * For @p{dim=3}. Adds (appends) the - * support points of all faces (quads in 3d) to - * the vector a. + * For @p{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 @p{dim=1} and 2 this - * function is empty. + * function. For @p{dim=1} and 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. */ //TODO: rename function to add_quad_support_points, to unify notation - void add_face_support_points(const typename Triangulation::cell_iterator &cell, - std::vector > &a) const; + virtual void + add_face_support_points(const typename Triangulation::cell_iterator &cell, + std::vector > &a) const; /** * For @p{dim=2} and 3. Simple @@ -375,6 +411,7 @@ class MappingQ : public MappingQ1 * Needed by the * @p{compute_support_points_simple} */ +//TODO: remove this function altogether? void fill_quad_support_points_simple (const Triangulation::cell_iterator &cell, std::vector > &a) const; @@ -445,7 +482,7 @@ class MappingQ : public MappingQ1 * Number of the Qp tensor * product shape functions. */ - unsigned int n_shape_functions; + const unsigned int n_shape_functions; /** * Mapping from lexicographic to diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index 304fe5c7a7..bd52d90ab5 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors +// Copyright (C) 2000, 2001 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -36,6 +36,7 @@ template class MappingQ1 : public Mapping { public: +//TODO: why make the following functions public? they are only helpful for fevalues and maybe the finite elements? /** * Implementation of the interface in * @ref{Mapping}. @@ -213,9 +214,10 @@ class MappingQ1 : public Mapping { public: /** - * Constructor. + * Constructor. Pass the + * number of shape functions. */ - InternalData(unsigned int n_shape_functions); + InternalData(const unsigned int n_shape_functions); /** * Shape function at quadrature @@ -224,31 +226,31 @@ class MappingQ1 : public Mapping * vertices must be reordered * to obtain transformation. */ - double shape (unsigned int qpoint, - unsigned int shape_nr) const; + double shape (const unsigned int qpoint, + const unsigned int shape_nr) const; /** * Shape function at quadrature * point. See above. */ - double &shape (unsigned int qpoint, - unsigned int shape_nr); + double &shape (const unsigned int qpoint, + const unsigned int shape_nr); /** * Gradient of shape function * in quadrature point. See * above. */ - Tensor<1,dim> derivative (unsigned int qpoint, - unsigned int shape_nr) const; + Tensor<1,dim> derivative (const unsigned int qpoint, + const unsigned int shape_nr) const; /** * Gradient of shape function * in quadrature point. See * above. */ - Tensor<1,dim> &derivative (unsigned int qpoint, - unsigned int shape_nr); + Tensor<1,dim> &derivative (const unsigned int qpoint, + const unsigned int shape_nr); /** * Values of shape @@ -299,11 +301,6 @@ class MappingQ1 : public Mapping */ std::vector > > aux; - /** - * Number of shape functions. - */ - unsigned int n_shape_functions; - /** * Stores the support points of * the mapping shape functions on @@ -316,7 +313,7 @@ class MappingQ1 : public Mapping * @p{mapping_support_points} are * stored. */ - DoFHandler::cell_iterator cell_of_current_support_points; + typename DoFHandler::cell_iterator cell_of_current_support_points; /** * Default value of this flag @@ -326,6 +323,21 @@ class MappingQ1 : public Mapping * @p{false}. */ bool is_mapping_q1_data; + + /** + * Number of shape + * functions. If this is a Q1 + * mapping, then it is simply + * the number of vertices per + * cell. However, since also + * derived classes use this + * class (e.g. the + * @ref{Mapping_Q} class), + * the number of shape + * functions may also be + * different. + */ + unsigned int n_shape_functions; }; /** @@ -424,9 +436,11 @@ class MappingQ1 : public Mapping std::vector > &a) const; /** - *Number of shape functions + * Number of shape functions. Is + * simply the number of vertices + * per cell for the Q1 mapping. */ - static const unsigned int n_shape_functions = 1 << dim; + static const unsigned int n_shape_functions = GeometryInfo::vertices_per_cell; }; diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index 269673debc..136e07ddcc 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -85,7 +85,7 @@ MappingQ::MappingQ (const unsigned int p): n_outer((dim==2) ? 4+4*(degree-1) :8+12*(degree-1)+6*(degree-1)*(degree-1)), tensor_pols(0), - n_shape_functions(0), + n_shape_functions(power(degree+1,dim)), renumber(0), //TODO: why have two ways to compute? if they both work, choose one and remove the other alternative_normals_computation(false), @@ -101,7 +101,8 @@ MappingQ::MappingQ (const unsigned int p): v.push_back(LagrangeEquidistant(degree,i)); tensor_pols = new TensorProductPolynomials (v); - n_shape_functions=tensor_pols->n_tensor_product_polynomials(); + Assert (n_shape_functions==tensor_pols->n_tensor_product_polynomials(), + ExcInternalError()); Assert(n_inner+n_outer==n_shape_functions, ExcInternalError()); // build the renumbering of the @@ -189,6 +190,7 @@ MappingQ::compute_shapes_virtual (const std::vector > &unit_poin } + template UpdateFlags MappingQ::update_each (const UpdateFlags in) const @@ -260,10 +262,10 @@ MappingQ::get_subface_data (const UpdateFlags update_flags, template void -MappingQ::compute_face_data (UpdateFlags update_flags, - const Quadrature& q, - const unsigned int n_original_q_points, - MappingQ1::InternalData& mapping_q1_data) const +MappingQ::compute_face_data (UpdateFlags update_flags, + const Quadrature &q, + const unsigned int n_original_q_points, + MappingQ1::InternalData &mapping_q1_data) const { // convert data object to internal // data for this class. fails with @@ -328,9 +330,17 @@ MappingQ::fill_fe_values (const DoFHandler::cell_iterator &cell, // possible InternalData &data = dynamic_cast (mapping_data); - data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells - || cell->has_boundary_lines()); - + // 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 if (data.use_mapping_q1_on_current_cell) MappingQ1::fill_fe_values(cell, q, data.mapping_q1_data, quadrature_points, JxW_values); @@ -340,6 +350,7 @@ MappingQ::fill_fe_values (const DoFHandler::cell_iterator &cell, } + template void MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, @@ -357,6 +368,11 @@ MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterato // possible InternalData &data = dynamic_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 //TODO: shouldn't we ask whether the face is at the boundary, rather than the cell? data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells || cell->has_boundary_lines()); @@ -364,6 +380,9 @@ MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterato const unsigned int npts=q.n_quadrature_points; const unsigned int offset=face_no*npts; + // depending on this result, use + // this or the other data object + // for the mapping if (data.use_mapping_q1_on_current_cell) MappingQ1::compute_fill_face (cell, face_no, false, npts, offset, q.get_weights(), @@ -397,6 +416,11 @@ MappingQ::fill_fe_subface_values (const typename DoFHandler::cell_iter // possible InternalData &data = dynamic_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 //TODO: shouldn't we ask whether the face is at the boundary, rather than the cell? data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells || cell->has_boundary_lines()); @@ -405,6 +429,9 @@ MappingQ::fill_fe_subface_values (const typename DoFHandler::cell_iter const unsigned int offset= (face_no*GeometryInfo::subfaces_per_face + sub_no)*npts; + // depending on this result, use + // this or the other data object + // for the mapping if (data.use_mapping_q1_on_current_cell) MappingQ1::compute_fill_face (cell, face_no, true, npts, offset, q.get_weights(), @@ -698,13 +725,18 @@ MappingQ::apply_laplace_vector(const std::vector > &lvs const unsigned int n_outer_apply=lvs[0].size(); 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 for (unsigned int unit_point=0; unit_point &lv=lvs[unit_point]; - Assert(lv.size()==n_outer_apply, ExcInternalError()); + Assert(lvs[unit_point].size()==n_outer_apply, ExcInternalError()); Point p; for (unsigned int k=0; k::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 (use_mapping_q_on_all_cells || cell->has_boundary_lines()) compute_support_points_laplace(cell, a); // compute_support_points_simple(cell, a); +//TODO: can we delete the previous line? + else + // 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); @@ -735,36 +778,43 @@ void MappingQ::compute_support_points_laplace(const typename Triangulation::cell_iterator &cell, std::vector > &a) const { + // in any case, we need the + // vertices first a.resize(GeometryInfo::vertices_per_cell); - // the vertices first for (unsigned int i=0; i::vertices_per_cell; ++i) a[i] = cell->vertex(i); if (degree>1) - { - if (dim==1) - { - Assert(false, ExcNotImplemented()); - } - else - { - // then the points on lines - // (for dim=2,3) - add_line_support_points (cell, a); - - if (dim==2) - apply_laplace_vector(laplace_on_quad_vector,a); - else if (dim==3) - { - add_face_support_points(cell, a); + switch (dim) + { + case 2: + // in 2d, add the + // points on the four + // bounding lines to + // the exterior (outer) + // points + add_line_support_points (cell, a); + apply_laplace_vector (laplace_on_quad_vector,a); + break; + + case 3: + // in 3d also add the + // points located on + // the boundary faces + add_line_support_points (cell, a); + add_face_support_points (cell, a); + apply_laplace_vector (laplace_on_hex_vector, a); + break; - apply_laplace_vector(laplace_on_hex_vector, a); - } - } - } + default 1: + Assert(false, ExcNotImplemented()); + break; + }; } + +//TODO: remove the following function altogether template void MappingQ::compute_support_points_simple(const typename Triangulation::cell_iterator &cell, @@ -869,38 +919,41 @@ void MappingQ::add_line_support_points (const Triangulation::cell_iterator &cell, std::vector > &a) const { - const Boundary *boundary = 0; - - std::vector > line_points; - if (degree>2) - line_points.resize(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 - for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) + // if we only need the midpoint, + // then ask for it. + if (degree==2) { - const typename Triangulation::line_iterator line = cell->line(line_no); - if (line->at_boundary()) - boundary=&line->get_triangulation().get_boundary(line->boundary_indicator()); - else - boundary=&straight_boundary; - - // if we only need the - // midpoint, then ask for - // it. otherwise call the more - // complicated functions - if (degree==2) + for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) a.push_back(boundary->get_new_point_on_line(line)); - else + } + else + // otherwise call the more + // complicated functions and ask + // for inner points from the + // boundary description + { + std::vector > 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 + for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) { + const typename Triangulation::line_iterator line = cell->line(line_no); + + const Boundary * const boundary + = (line->at_boundary() ? + &line->get_triangulation().get_boundary(line->boundary_indicator()) : + &straight_boundary); + boundary->get_intermediate_points_on_line (line, line_points); a.insert (a.end(), line_points.begin(), line_points.end()); - } + } } -} +}; @@ -911,7 +964,7 @@ MappingQ::add_line_support_points (const Triangulation::cell_iterator template<> void MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell, - std::vector > &a) const + std::vector > &a) const { const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell, vertices_per_face = GeometryInfo<3>::vertices_per_face, @@ -956,6 +1009,7 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell // on it if (face->at_boundary()) { +//TODO: move this variable out of the inner loop std::vector > quad_points ((degree-1)*(degree-1)); face->get_triangulation().get_boundary(face->boundary_indicator()) @@ -1011,6 +1065,7 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell // intermediate points // from a straight // boundary object +//TODO: move this variable out of the loop std::vector > quad_points ((degree-1)*(degree-1)); straight_boundary.get_intermediate_points_on_quad (face, quad_points); diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index 29713e2389..19a896fa1c 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -26,17 +26,22 @@ +template +const unsigned int MappingQ1::n_shape_functions; + + + template -MappingQ1::InternalData::InternalData (unsigned int n_shape_functions): - n_shape_functions(n_shape_functions), - is_mapping_q1_data(true) +MappingQ1::InternalData::InternalData (const unsigned int n_shape_functions): + is_mapping_q1_data(true), + n_shape_functions (n_shape_functions) {} template inline double -MappingQ1::InternalData::shape (unsigned int qpoint, - unsigned int shape_nr) const +MappingQ1::InternalData::shape (const unsigned int qpoint, + const unsigned int shape_nr) const { Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(), ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, shape_values.size())); @@ -300,10 +305,10 @@ MappingQ1::update_each (const UpdateFlags in) const template void -MappingQ1::compute_data (const UpdateFlags update_flags, - const Quadrature& q, - const unsigned int n_original_q_points, - InternalData& data) const +MappingQ1::compute_data (const UpdateFlags update_flags, + const Quadrature &q, + const unsigned int n_original_q_points, + InternalData &data) const { const unsigned int npts = q.n_quadrature_points; @@ -313,8 +318,6 @@ MappingQ1::compute_data (const UpdateFlags update_flags, const UpdateFlags flags(data.update_flags); - // cerr << "Data: " << hex << flags << dec << endl; - if (flags & update_transformation_values) data.shape_values.resize(data.n_shape_functions * npts); @@ -411,8 +414,8 @@ MappingQ1::compute_face_data (UpdateFlags update_flags, template Mapping::InternalDataBase* -MappingQ1::get_face_data (const UpdateFlags update_flags, - const Quadrature& quadrature) const +MappingQ1::get_face_data (const UpdateFlags update_flags, + const Quadrature &quadrature) const { InternalData* data = new InternalData(n_shape_functions); QProjector q (quadrature, false); @@ -1022,7 +1025,7 @@ MappingQ1::covariant_transformation (std::vector &dst, typename std::vector::const_iterator vec = src.begin() + src_offset; typename std::vector >::const_iterator tensor = data.covariant.begin(); typename std::vector::iterator result = dst.begin(); - typename std::vector::const_iterator end = dst.end(); + const typename std::vector::const_iterator end = dst.end(); while (result!=end) { -- 2.39.5