From: Wolfgang Bangerth Date: Thu, 8 Mar 2001 10:35:32 +0000 (+0000) Subject: Clean-ups, annotations, etc. X-Git-Tag: v8.0.0~19601 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9909ed0ec95fd7265bb36adb463dfdb7d8c7ee62;p=dealii.git Clean-ups, annotations, etc. git-svn-id: https://svn.dealii.org/trunk@4170 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index fed9d47ca0..dfc0513c51 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -111,6 +111,7 @@ class Mapping : public Subscriptor * Virtual destructor. */ virtual ~Mapping (); + /** * Prepare internal data * structures and fill in values @@ -129,6 +130,7 @@ class Mapping : public Subscriptor virtual InternalDataBase* get_face_data (const UpdateFlags flags, const Quadrature& quadrature) const = 0; + /** * Prepare internal data * structure for transformation diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h index be6c2e2374..5b5d82bc13 100644 --- a/deal.II/deal.II/include/fe/mapping_cartesian.h +++ b/deal.II/deal.II/include/fe/mapping_cartesian.h @@ -142,6 +142,7 @@ class MappingCartesian : public Mapping * @p{p_real} on the real cell * @p{cell} and returns @p{p_real}. */ +//TODO: document meaning of mdata argument virtual Point transform_unit_to_real_cell ( const typename Triangulation::cell_iterator cell, const Point &p, @@ -209,14 +210,14 @@ class MappingCartesian : public Mapping /** * Constructor. */ - InternalData (const Quadrature& quadrature); + InternalData (const Quadrature &quadrature); /** * Length of the cell in * different coordinate * directions. */ - std::vector length; + Tensor<1,dim> length; /** * Vector of all quadrature @@ -236,7 +237,7 @@ class MappingCartesian : public Mapping std::vector > > unit_tangentials; /** - * Auxuliary vectors for internal use. + * Auxiliary vectors for internal use. */ std::vector > > aux; }; @@ -252,6 +253,14 @@ class MappingCartesian : public Mapping InternalData& data, std::vector > &quadrature_points, std::vector >& normal_vectors) const; + + private: + /** + * Value to indicate that a given + * face or subface number is + * invalid. + */ + static const unsigned int invalid_face_number = static_cast(-1); }; diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index 5c33194d56..4938a56974 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -25,7 +25,11 @@ class LagrangeEquidistant; //TODO: fill_fe_face_values should exist in a version doing all faces -// to save initialization time. +//TODO: to save initialization time. + +//TODO: doc: laplace_on_quad_vector, compute_laplace_on_quad, etc all +//TODO: reference to each other, there is no description what they actually +//TODO: do or are good for /** * Mapping class that uses Qp-mappings on boundary AND on inner @@ -269,7 +273,7 @@ class MappingQ : public MappingQ1 * * For @p{degree==2} this function * sets the - * @p{laplace_on_quad_vector} to + * @p{laplace_on_hex_vector} to * the hardcoded data. For * @p{degree>2} this vector is * computed. @@ -358,6 +362,7 @@ class MappingQ : public MappingQ1 * function. For @p{dim=1} and 2 this * function is empty. */ +//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; diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index 46450cfae4..304fe5c7a7 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -147,6 +147,7 @@ class MappingQ1 : public Mapping * @p{p_real} on the real cell * @p{cell} and returns @p{p_real}. */ +//TODO: document meaning of mdata argument virtual Point transform_unit_to_real_cell ( const typename Triangulation::cell_iterator cell, const Point &p, diff --git a/deal.II/deal.II/include/grid/tria_boundary.h b/deal.II/deal.II/include/grid/tria_boundary.h index 4503cd8030..80ab238564 100644 --- a/deal.II/deal.II/include/grid/tria_boundary.h +++ b/deal.II/deal.II/include/grid/tria_boundary.h @@ -121,9 +121,11 @@ class Boundary : public Subscriptor get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; /** - * Return intermediate points - * on the boundary line. - * + * Return equally spaced + * intermediate points on a + * boundary line. +//TODO: this function is also called for inner lines!? + * * The number of points requested * is given by the size of the * vector @p{points}. It is the @@ -145,8 +147,10 @@ class Boundary : public Subscriptor typename std::vector > &points) const; /** - * Return intermediate points - * on the boundary quad. + * Return equally spaced + * intermediate points on a + * boundary quad. +//TODO: this function is also called for inner quads!? * * The number of points requested * is given by the size of the @@ -235,8 +239,9 @@ class StraightBoundary : public Boundary { * p{n+1} partitions of equal * lengths. * - * Refer to the general documentation of - * this class and the documentation of the + * Refer to the general + * documentation of this class + * and the documentation of the * base class. */ virtual void @@ -250,8 +255,9 @@ class StraightBoundary : public Boundary { * @p{(m+1)(m+1)} subquads of equal * size. * - * Refer to the general documentation of - * this class and the documentation of the + * Refer to the general + * documentation of this class + * and the documentation of the * base class. */ virtual void diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index 670f8a122d..e9d98c7a3f 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -60,7 +60,9 @@ FE_DGQ::FE_DGQ (unsigned int degree) poly = new TensorProductPolynomials (v); } - // generate permutation/rotation index sets to generate some matrices from others + // generate permutation/rotation + // index sets to generate some + // matrices from others std::vector right; std::vector top; rotate_indices (right, 'Z'); diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 449145ec4a..90b17638d2 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -49,6 +49,7 @@ FE_Q::FE_Q (const unsigned int degree) // do some internal book-keeping build_renumbering (*this, degree, renumber); +//TODO: externalize this to a proper template function #if deal_II_dimension > 1 build_face_renumbering (FiniteElementData(FE_Q::get_dpo_vector(degree),1), degree, diff --git a/deal.II/deal.II/source/fe/mapping.cc b/deal.II/deal.II/source/fe/mapping.cc index 8452ebf17e..0a8fd3918d 100644 --- a/deal.II/deal.II/source/fe/mapping.cc +++ b/deal.II/deal.II/source/fe/mapping.cc @@ -18,7 +18,7 @@ #if (deal_II_dimension == 1) template<> -const unsigned int Mapping::normal_directions[2] = +const unsigned int Mapping<1>::normal_directions[2] = { 1, 0 }; @@ -28,7 +28,7 @@ const unsigned int Mapping::normal_directions[2] = #if (deal_II_dimension == 2) template<> -const unsigned int Mapping::normal_directions[4] = +const unsigned int Mapping<2>::normal_directions[4] = { 2, 0, 3, 1 }; @@ -38,7 +38,7 @@ const unsigned int Mapping::normal_directions[4] = #if (deal_II_dimension == 3) template<> -const unsigned int Mapping::normal_directions[6] = +const unsigned int Mapping<3>::normal_directions[6] = { 3, 2, 5, 0, 4, 1 }; @@ -51,12 +51,7 @@ Mapping::~Mapping () {} -/*------------------------------ InternalData ------------------------------*/ - - -template -Mapping::InternalDataBase::~InternalDataBase () -{} +/*------------------------------ InternalDataBase ------------------------------*/ template @@ -69,6 +64,11 @@ Mapping::InternalDataBase::InternalDataBase (): +template +Mapping::InternalDataBase::~InternalDataBase () +{} + + /*------------------------------ InternalData ------------------------------*/ diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index 49fad8d038..805a8a4515 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -26,10 +26,14 @@ +template +const unsigned int MappingCartesian::invalid_face_number; + + + template MappingCartesian::InternalData::InternalData (const Quadrature& q) : - length (dim, 0.), quadrature_points (q.get_points ()) {} @@ -42,6 +46,8 @@ MappingCartesian::update_once (const UpdateFlags) const return update_default; } + + template UpdateFlags MappingCartesian::update_each (const UpdateFlags in) const @@ -54,10 +60,11 @@ MappingCartesian::update_each (const UpdateFlags in) const } + template Mapping::InternalDataBase* -MappingCartesian::get_data (const UpdateFlags flags, - const Quadrature& q) const +MappingCartesian::get_data (const UpdateFlags flags, + const Quadrature &q) const { Assert (flags & update_normal_vectors == 0, ExcNotImplemented()); InternalData* data = new InternalData (q); @@ -68,7 +75,7 @@ MappingCartesian::get_data (const UpdateFlags flags, template Mapping::InternalDataBase* MappingCartesian::get_face_data (const UpdateFlags, - const Quadrature& quadrature) const + const Quadrature& quadrature) const { QProjector q (quadrature, false); InternalData* data = new InternalData (q); @@ -80,7 +87,7 @@ MappingCartesian::get_face_data (const UpdateFlags, template Mapping::InternalDataBase* MappingCartesian::get_subface_data (const UpdateFlags, - const Quadrature& quadrature) const + const Quadrature &quadrature) const { QProjector q (quadrature, true); InternalData* data = new InternalData (q); @@ -93,45 +100,55 @@ MappingCartesian::get_subface_data (const UpdateFlags, template void MappingCartesian::compute_fill (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - InternalData& data, + const unsigned int face_no, + const unsigned int sub_no, + InternalData &data, std::vector > &quadrature_points, - std::vector >& normal_vectors) const + std::vector > &normal_vectors) const { - const UpdateFlags update_flags(data.current_update_flags()); +//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions? + UpdateFlags update_flags(data.current_update_flags()); const unsigned int npts = quadrature_points.size (); unsigned int offset = 0; - bool onface = false; - if (face_no == static_cast (-1)) + if (face_no != invalid_face_number) { - Assert (sub_no == static_cast (-1), ExcInternalError()); - } else { - onface = true; // Add 1 on both sides of // assertion to avoid compiler // warning about testing // unsigned int < 0 in 1d. Assert (face_no+1 < GeometryInfo::faces_per_cell+1, ExcIndexRange (face_no, 0, GeometryInfo::faces_per_cell)); - if (sub_no == static_cast (-1)) + + if (sub_no == invalid_face_number) + // called from FEFaceValues offset = face_no * quadrature_points.size(); else { + // called from FESubfaceValue Assert (sub_no+1 < GeometryInfo::subfaces_per_face+1, ExcIndexRange (sub_no, 0, GeometryInfo::subfaces_per_face)); offset = (face_no * GeometryInfo::subfaces_per_face + sub_no) * quadrature_points.size(); } } + else + // invalid face number, so + // subface should be invalid as + // well + Assert (sub_no == invalid_face_number, ExcInternalError()); + + // let @p{start} be the origin of a + // local coordinate system. it is + // chosen as the (lower) left + // vertex + const Point start = cell->vertex(0); - // Compute start point and sizes along axes. - // Strange vertex numbering makes this complicated again. - - Point start = cell->vertex (0); - + // Compute start point and sizes + // along axes. Strange vertex + // numbering makes this complicated + // again. switch (dim) { case 1: @@ -150,23 +167,34 @@ MappingCartesian::compute_fill (const typename DoFHandler::cell_iterat Assert(false, ExcNotImplemented()); } - + + // transform quadrature point. this + // is obtained simply by scaling + // unit coordinates with lengths in + // each direction if (update_flags & update_q_points) { Assert (quadrature_points.size() == npts, ExcDimensionMismatch(quadrature_points.size(), npts)); - for (unsigned int i=0;i p = start; - for (unsigned int d=0;d n; switch (100*dim+face_no) { @@ -205,8 +233,9 @@ MappingCartesian::compute_fill (const typename DoFHandler::cell_iterat default: Assert (false, ExcInternalError()); } - for (unsigned int i=0;i::fill_fe_values (const DoFHandler::cell_iterator& cel std::vector >& quadrature_points, std::vector& JxW_values) const { - InternalData *data_ptr = dynamic_cast (&mapping_data); - Assert(data_ptr!=0, ExcInternalError()); - InternalData &data=*data_ptr; + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &data = dynamic_cast (mapping_data); typename std::vector > dummy; - compute_fill (cell, static_cast (-1), static_cast (-1), + compute_fill (cell, invalid_face_number, invalid_face_number, data, quadrature_points, dummy); + // compute Jacobian + // determinant. all values are + // equal and are the product of the + // local lengths in each coordinate + // direction +//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions? if (data.current_update_flags() & update_JxW_values) { double J = data.length[0]; @@ -247,34 +284,42 @@ MappingCartesian::fill_fe_values (const DoFHandler::cell_iterator& cel template void MappingCartesian::fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const Quadrature &q, - typename Mapping::InternalDataBase &mapping_data, - std::vector > &quadrature_points, - std::vector &JxW_values, - std::vector > &boundary_forms, - std::vector > &normal_vectors) const + const unsigned int face_no, + const Quadrature &q, + typename Mapping::InternalDataBase &mapping_data, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &boundary_forms, + std::vector > &normal_vectors) const { - InternalData *data_ptr = dynamic_cast (&mapping_data); - Assert(data_ptr!=0, ExcInternalError()); - InternalData &data=*data_ptr; + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &data = dynamic_cast (mapping_data); - compute_fill (cell, face_no, static_cast (-1), + compute_fill (cell, face_no, invalid_face_number, data, quadrature_points, normal_vectors); + // first compute Jacobian + // determinant, which is simply the + // product of the local lengths + // since the jacobian is diagonal double J = 1.; for (unsigned int d=0;d::fill_fe_subface_values (const typename DoFHandler::c std::vector > &boundary_forms, std::vector > &normal_vectors) const { - InternalData *data_ptr = dynamic_cast (&mapping_data); - Assert(data_ptr!=0, ExcInternalError()); - InternalData &data=*data_ptr; + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &data = dynamic_cast (mapping_data); compute_fill (cell, face_no, sub_no, data, quadrature_points, normal_vectors); + // first compute Jacobian + // determinant, which is simply the + // product of the local lengths + // since the jacobian is diagonal double J = 1.; for (unsigned int d=0;d::subfaces_per_face; } +//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions? if (data.current_update_flags() & update_boundary_forms) { for (unsigned int i=0; i::fill_fe_subface_values (const typename DoFHandler::c template <> void MappingCartesian<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &, - const unsigned, - const Quadrature<0>&, - Mapping<1>::InternalDataBase&, - std::vector >&, - std::vector&, - std::vector >&, - std::vector >&) const + const unsigned, + const Quadrature<0>&, + Mapping<1>::InternalDataBase&, + std::vector >&, + std::vector&, + std::vector >&, + std::vector >&) const { Assert(false, ExcNotImplemented()); } @@ -343,14 +396,14 @@ MappingCartesian<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &, template <> void MappingCartesian<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterator &, - const unsigned, - const unsigned, - const Quadrature<0>&, - Mapping<1>::InternalDataBase&, - std::vector >&, - std::vector&, - std::vector >&, - std::vector >&) const + const unsigned, + const unsigned, + const Quadrature<0>&, + Mapping<1>::InternalDataBase&, + std::vector >&, + std::vector&, + std::vector >&, + std::vector >&) const { Assert(false, ExcNotImplemented()); } @@ -407,6 +460,8 @@ Point MappingCartesian::transform_unit_to_real_cell ( const Point &p, const typename Mapping::InternalDataBase *const m_data) const { +//TODO: wouldn't it be simpler in this function to compute local lengths ourselved, rather than relying on an InternalDataBase object? + // If m_data!=0 use this // InternalData. // @@ -432,12 +487,19 @@ Point MappingCartesian::transform_unit_to_real_cell ( mdata = dynamic_cast (m_data); Assert(mdata!=0, ExcInternalError()); +//TODO: don't we need to reinit() mdata on the present cell or initialize the mdata->length array in some other way? + // use now the InternalData, that // mdata is pointing to, to compute - // the point in real space. + // the point in real space. it is + // simply a scaled version of the + // unit cell point shifted by the + // (lower) left corner of the cell Point p_real = cell->vertex(0); - for (unsigned int d=0;dlength[d]*p(d); + +//TODO: memory leak: if initially mdata==0, then an object is created but never deleted. see mapping_q1 for a better implementation return p_real; } @@ -448,7 +510,7 @@ Point MappingCartesian::transform_real_to_unit_cell ( const typename Triangulation::cell_iterator cell, const Point &p) const { - const Point& start = cell->vertex (0); + const Point &start = cell->vertex(0); Point real = p; real -= start; @@ -479,53 +541,64 @@ template inline void MappingCartesian::contravariant_transformation (std::vector &dst, - const std::vector &src, - const Mapping::InternalDataBase &mapping_data, - const unsigned int src_offset) const + const std::vector &src, + const Mapping::InternalDataBase &mapping_data, + const unsigned int src_offset) const { Assert(tensor_::dimension==dim && tensor_::rank==1, ExcInvalidData()); - const InternalData* data_ptr = dynamic_cast (&mapping_data); - Assert(data_ptr!=0, ExcInternalError()); - const InternalData &data=*data_ptr; + + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + const InternalData &data = dynamic_cast (mapping_data); Assert (data.update_flags & update_contravariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - typename std::vector::const_iterator vec = src.begin() + src_offset; - typename std::vector::iterator result = dst.begin(); - typename std::vector::const_iterator end = dst.end(); + typename std::vector::const_iterator vec = src.begin() + src_offset; + typename std::vector::iterator result = dst.begin(); + typename std::vector::const_iterator end = dst.end(); + // simply scale by Jacobian + // (which is diagonal here) while (result!=end) { - for (unsigned int d=0;d template inline void MappingCartesian::covariant_transformation (std::vector &dst, - const std::vector &src, - const Mapping::InternalDataBase &mapping_data, - const unsigned int src_offset) const + const std::vector &src, + const Mapping::InternalDataBase &mapping_data, + const unsigned int src_offset) const { Assert(tensor_::dimension==dim && tensor_::rank==1, ExcInvalidData()); - const InternalData *data_ptr = dynamic_cast (&mapping_data); - Assert(data_ptr!=0, ExcInternalError()); - const InternalData &data=*data_ptr; + + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + const InternalData &data = dynamic_cast (mapping_data); Assert (data.update_flags & update_covariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - typename std::vector::const_iterator vec = src.begin() + src_offset; - typename std::vector::iterator result = dst.begin(); - typename std::vector::const_iterator end = dst.end(); - + typename std::vector::const_iterator vec = src.begin() + src_offset; + typename std::vector::iterator result = dst.begin(); + typename std::vector::const_iterator end = dst.end(); + + // simply scale by inverse Jacobian + // (which is diagonal here) while (result!=end) { for (unsigned int d=0;d::covariant_transformation (std::vector &dst //----------------------------------------------------------------------// +// explicit instantiations template class MappingCartesian; diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index f42fe23efb..269673debc 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -27,10 +27,11 @@ template -MappingQ::InternalData::InternalData (unsigned int n_shape_functions): +MappingQ::InternalData::InternalData (const unsigned int n_shape_functions): MappingQ1::InternalData(n_shape_functions), - use_mapping_q1_on_current_cell(false), - mapping_q1_data(1 << dim) +//TODO: in 1d, use_mapping_q1_on_current_cell is always false, but should be true. + use_mapping_q1_on_current_cell(false), + mapping_q1_data(1 << dim) { is_mapping_q1_data=false; } @@ -38,8 +39,10 @@ MappingQ::InternalData::InternalData (unsigned int n_shape_functions): #if deal_II_dimension == 1 +// in 1d, it is irrelevant which polynomial degree to use, since all +// cells are scaled linearly template<> -MappingQ<1>::MappingQ (unsigned int): +MappingQ<1>::MappingQ (const unsigned int): laplace_on_quad_vector(0), laplace_on_hex_vector(0), degree(1), @@ -94,7 +97,7 @@ MappingQ::MappingQ (const unsigned int p): // functions for the Qp mapping of // cells at the boundary. std::vector v; - for (unsigned int i=0;i<=degree;++i) + for (unsigned int i=0; i<=degree; ++i) v.push_back(LagrangeEquidistant(degree,i)); tensor_pols = new TensorProductPolynomials (v); @@ -107,15 +110,16 @@ MappingQ::MappingQ (const unsigned int p): renumber.resize(n_shape_functions,0); std::vector dpo(dim+1, 1); for (unsigned int i=1; i fe_data(dpo, 1); FE_Q::build_renumbering (fe_data, p, renumber); // build laplace_on_quad_vector if (degree>1) { - set_laplace_on_quad_vector(laplace_on_quad_vector); - if (dim==3) + if (dim >= 2) + set_laplace_on_quad_vector(laplace_on_quad_vector); + if (dim >= 3) set_laplace_on_hex_vector(laplace_on_hex_vector); } } @@ -134,7 +138,7 @@ MappingQ::~MappingQ () template<> void MappingQ<1>::compute_shapes_virtual (const std::vector > &unit_points, - MappingQ1<1>::InternalData &data) const + MappingQ1<1>::InternalData &data) const { MappingQ1<1>::compute_shapes_virtual(unit_points, data); } @@ -199,6 +203,7 @@ MappingQ::update_each (const UpdateFlags in) const } + template Mapping::InternalDataBase* MappingQ::get_data (const UpdateFlags update_flags, @@ -260,13 +265,16 @@ MappingQ::compute_face_data (UpdateFlags update_flags, const unsigned int n_original_q_points, MappingQ1::InternalData& mapping_q1_data) const { - InternalData *data_ptr = dynamic_cast (&mapping_q1_data); - Assert(data_ptr!=0, ExcInternalError()); - InternalData &data=*data_ptr; + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &data = dynamic_cast (mapping_q1_data); MappingQ1::compute_face_data(update_flags, q, n_original_q_points, data); - + +//TODO: externalize this to a proper template function. or: scrap the whole function for 1d since it doesn't make much sense anyway? #if (deal_II_dimension>1) if ((data.update_flags & update_normal_vectors) && alternative_normals_computation) @@ -310,13 +318,15 @@ template void MappingQ::fill_fe_values (const DoFHandler::cell_iterator &cell, const Quadrature &q, - Mapping::InternalDataBase &mapping_data, - std::vector > &quadrature_points, - std::vector &JxW_values) const + Mapping::InternalDataBase &mapping_data, + std::vector > &quadrature_points, + std::vector &JxW_values) const { - InternalData *data_ptr = dynamic_cast (&mapping_data); - Assert(data_ptr!=0, ExcInternalError()); - InternalData &data=*data_ptr; + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &data = dynamic_cast (mapping_data); data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells || cell->has_boundary_lines()); @@ -341,10 +351,13 @@ MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterato std::vector > &exterior_forms, std::vector > &normal_vectors) const { - InternalData *data_ptr = dynamic_cast (&mapping_data); - Assert(data_ptr!=0, ExcInternalError()); - InternalData &data=*data_ptr; + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &data = dynamic_cast (mapping_data); +//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()); @@ -378,10 +391,13 @@ MappingQ::fill_fe_subface_values (const typename DoFHandler::cell_iter std::vector > &exterior_forms, std::vector > &normal_vectors) const { - InternalData *data_ptr = dynamic_cast (&mapping_data); - Assert(data_ptr!=0, ExcInternalError()); - InternalData &data=*data_ptr; + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &data = dynamic_cast (mapping_data); +//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()); @@ -422,7 +438,10 @@ MappingQ::set_laplace_on_quad_vector(std::vector > &loq Assert(degree>1, ExcInternalError()); 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 double const *loqv_ptr=0; if (degree==2) { @@ -432,7 +451,7 @@ MappingQ::set_laplace_on_quad_vector(std::vector > &loq } else if (degree==3) { - static double loqv3[4*12] + static const double loqv3[4*12] ={80/1053., 1/81., 11/1053., 1/81., 25/117., 44/351., 7/117., 16/351., 7/117., 16/351., 25/117., 44/351., 1/81., 80/1053., 1/81., 11/1053., 44/351., 25/117., @@ -444,29 +463,33 @@ MappingQ::set_laplace_on_quad_vector(std::vector > &loq loqv_ptr=&loqv3[0]; } - + if (loqv_ptr!=0) { + // precomputed. copy values to + // the loqvs array loqvs.resize(n_inner_2d); for (unsigned int unit_point=0; unit_point &loqv=loqvs[unit_point]; for (unsigned int k=0; k::set_laplace_on_hex_vector(std::vector > &lohvs) { Assert(degree>1, ExcInternalError()); + // first check whether we have + // precomputed the values for some + // polynomial degree double const *lohv_ptr=0; if (degree==2) { @@ -504,22 +530,26 @@ MappingQ<3>::set_laplace_on_hex_vector(std::vector > &lohvs) if (lohv_ptr!=0) { + // precomputed. copy values to + // the lohvs array lohvs.resize(n_inner); for (unsigned int unit_point=0; unit_point &loqv=lohvs[unit_point]; for (unsigned int k=0; k::compute_laplace_vector(std::vector > &lvs) co // compute the shape // gradients at the quadrature // points on the unit cell - QGauss4 quadrature; + const QGauss4 quadrature; const unsigned int n_q_points=quadrature.n_quadrature_points; InternalData quadrature_data(n_shape_functions); - quadrature_data.shape_derivatives - .resize(n_shape_functions * n_q_points); + quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points); compute_shapes(quadrature.get_points(), quadrature_data); // Compute the stiffness matrix of @@ -706,10 +735,10 @@ void MappingQ::compute_support_points_laplace(const typename Triangulation::cell_iterator &cell, std::vector > &a) const { - a.resize(0); + a.resize(GeometryInfo::vertices_per_cell); // the vertices first for (unsigned int i=0; i::vertices_per_cell; ++i) - a.push_back(cell->vertex(i)); + a[i] = cell->vertex(i); if (degree>1) { @@ -742,6 +771,7 @@ MappingQ::compute_support_points_simple(const typename Triangulation:: std::vector > &a) const { // the vertices first +// TODO: check for size of 'a' first? for (unsigned int i=0; i::vertices_per_cell; ++i) a.push_back(cell->vertex(i)); @@ -755,64 +785,67 @@ MappingQ::compute_support_points_simple(const typename Triangulation:: // (for dim=3) fill_quad_support_points_simple (cell, a); - // then the points in cell - Point middle; - compute_midpoint(a, middle); - - if (degree==2) - a.push_back(middle); - else if (degree==3) - // The four points in the - // cell are located at the - // midpoint between the - // middle point and the 4 - // vertices - - // TODO: better position of - // points: transform them by - // a Q2 transformation. - for (unsigned int i=0; i::vertices_per_cell; ++i) - a.push_back(middle*2./3.+cell->vertex(vertex_mapping[i])/3.); - else if (degree==4) - { - Assert(a.size()==16, ExcInternalError()); - a.insert(a.end(), 9, Point()); + // then the points in cell. for + // this we need the midpoint of + // the points already in @p{a} + const Point middle = std::accumulate(a.begin(), a.end(), Point()) + / a.size(); - const unsigned int inner_map[8]= - { 0, 1, 2, 5, 8, 7, 6, 3 }; - - - // The nine points in the - // cell are located at the - // midpoint between the - // middle point and (the 4 - // vertices and the face - // midpoints) - - a[16+4]=middle; - for (unsigned int i=0, j=0; i::vertices_per_cell; ++i) - { - a[16+inner_map[j++]]=(middle+cell->vertex(i))/2.; - a[16+inner_map[j++]]=(middle+(cell->vertex(i)+cell->vertex((i+1)%4))/2.)/2.; - } - } - else - Assert(false, ExcNotImplemented()); - } + switch (degree) + { + case 2: + { + a.push_back(middle); + break; + }; + + case 3: + { + // The four points in the + // cell are located at + // the midpoint between + // the middle point and + // the 4 vertices + + // TODO: better position of + // points: transform them by + // a Q2 transformation. + for (unsigned int i=0; i::vertices_per_cell; ++i) + a.push_back(middle*2./3.+cell->vertex(vertex_mapping[i])/3.); + break; + }; + + case 4: + { + Assert(a.size()==16, ExcInternalError()); + a.insert(a.end(), 9, Point()); + + const unsigned int inner_map[8]= + { 0, 1, 2, 5, 8, 7, 6, 3 }; + + + // The nine points in the + // cell are located at the + // midpoint between the + // middle point and (the 4 + // vertices and the face + // midpoints) + + a[16+4]=middle; + for (unsigned int i=0, j=0; i::vertices_per_cell; ++i) + { + a[16+inner_map[j++]]=(middle+cell->vertex(i))/2.; + a[16+inner_map[j++]]=(middle+(cell->vertex(i)+cell->vertex((i+1)%4))/2.)/2.; + } + break; + }; + + default: + Assert(false, ExcNotImplemented()); + }; + }; + Assert(a.size()==n_shape_functions, ExcInternalError()); - -} - - - -template static -void compute_midpoint(const std::vector > &points, Point &m) -{ - typename std::vector >::const_iterator p_iter=points.begin(); - const typename std::vector >::const_iterator p_end=points.end(); - for (; p_iter!=p_end; ++p_iter) - m+=*p_iter; - m/=points.size(); } @@ -822,8 +855,11 @@ void compute_midpoint(const std::vector > &points, Point &m) template <> void MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &, - std::vector > &) const -{} + std::vector > &) const +{ + // there are no points on bounding + // lines which are to be added +} #endif @@ -831,30 +867,37 @@ MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &, template void MappingQ::add_line_support_points (const Triangulation::cell_iterator &cell, - std::vector > &a) const + std::vector > &a) const { - const Boundary *boundary; + const Boundary *boundary = 0; std::vector > line_points; if (degree>2) line_points.resize(degree-1); - - Triangulation::line_iterator line; + + // 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) { - line = cell->line(line_no); + 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) a.push_back(boundary->get_new_point_on_line(line)); else { boundary->get_intermediate_points_on_line (line, line_points); - for (unsigned int i=0; i::add_line_support_points (const Triangulation::cell_iterator #if deal_II_dimension==3 +//TODO: rename function to add_quad_support_points, to unify notation 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, - lines_per_face=GeometryInfo<3>::lines_per_face, - vertices_per_cell=GeometryInfo<3>::vertices_per_cell; + const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell, + vertices_per_face = GeometryInfo<3>::vertices_per_face, + lines_per_face = GeometryInfo<3>::lines_per_face, + vertices_per_cell = GeometryInfo<3>::vertices_per_cell; static const unsigned int face_vertex_to_cell_vertex [faces_per_cell][vertices_per_face]={{0,1,2,3}, @@ -891,10 +935,11 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell {8,7,11,3}}; - Triangulation<3>::face_iterator face; + // loop over all faces and collect points on them for (unsigned int face_no=0; face_noface(face_no); + const Triangulation<3>::face_iterator face=cell->face(face_no); + for (unsigned int i=0; ivertex_index(i)== cell->vertex_index(face_vertex_to_cell_vertex[face_no][i]), @@ -904,30 +949,41 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell Assert(face->line(i)== cell->line(face_line_to_cell_line[face_no][i]), ExcInternalError()); - + + // if face at boundary, then + // ask boundary object to + // return intermediate points + // on it if (face->at_boundary()) { - std::vector > quad_points; - quad_points.resize((degree-1)*(degree-1)); + std::vector > quad_points ((degree-1)*(degree-1)); face->get_triangulation().get_boundary(face->boundary_indicator()) .get_intermediate_points_on_quad (face, quad_points); - - for (unsigned int i=0; iline(i)->at_boundary()) ++lines_at_boundary; Assert(lines_at_boundary0) { // sort the points into b +//TODO: this is not thread-safe!!! b might be used for objects with +//TODO: different degrees at the same time! static std::vector > b; b.resize(4*degree); Assert(4*degree==vertices_per_face+lines_per_face*(degree-1), @@ -950,13 +1006,15 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell } else { - std::vector > quad_points; - quad_points.resize((degree-1)*(degree-1)); + // face is entirely in + // the interior. get + // intermediate points + // from a straight + // boundary object + std::vector > quad_points ((degree-1)*(degree-1)); straight_boundary.get_intermediate_points_on_quad (face, quad_points); - - for (unsigned int i=0; i::add_face_support_points(const typename Triangulation::cell_i template <> void MappingQ<3>::fill_quad_support_points_simple (const Triangulation<3>::cell_iterator &cell, - std::vector > &a) const + std::vector > &a) const { - const Boundary<3> *boundary; + const Boundary<3> *boundary = 0; std::vector > quad_points; Assert(degree>1, ExcInternalError()); quad_points.resize((degree-1)*(degree-1)); - Triangulation<3>::quad_iterator quad; for (unsigned int quad_no=0; quad_no::quads_per_cell; ++quad_no) { - quad = cell->face(quad_no); + const Triangulation<3>::quad_iterator quad = cell->face(quad_no); if (quad->at_boundary()) boundary=&quad->get_triangulation().get_boundary(quad->boundary_indicator()); else boundary=&straight_boundary; boundary->get_intermediate_points_on_quad (quad, quad_points); - for (unsigned int i=0; i::transform_contravariant (std::vector > &dst, - +// explicit instantiation template class MappingQ; diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index 79b70e7c79..29713e2389 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -864,6 +864,7 @@ Point MappingQ1::transform_unit_to_real_cell ( mdata=mdata_local; } else +//TODO: can we assume that m_data->mapping_support_points is initialized mdata = dynamic_cast (m_data); Assert(mdata!=0, ExcInternalError());