From: Wolfgang Bangerth Date: Wed, 22 Feb 2006 07:00:03 +0000 (+0000) Subject: Completely rewrite the base class for the hp::FE*Values classes to allow for arbitrar... X-Git-Tag: v8.0.0~12227 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7b78f7aa2c55f55eaa95a2ebb05dcc7282ecdba7;p=dealii.git Completely rewrite the base class for the hp::FE*Values classes to allow for arbitrary combinations of finite elements, quadrature objects, and mappings. git-svn-id: https://svn.dealii.org/trunk@12455 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/hp_fe_values.h b/deal.II/deal.II/include/fe/hp_fe_values.h index e2d6871c39..19d86c94cf 100644 --- a/deal.II/deal.II/include/fe/hp_fe_values.h +++ b/deal.II/deal.II/include/fe/hp_fe_values.h @@ -32,130 +32,23 @@ namespace internal namespace hp { /** - * Map between finite element objects and @p FEValues (the second - * template parameter, which might be FEValues, - * FEFaceValues, ...) objects. The hpFE*Values classes - * use this to hold an FE*Values object for each finite element - * that is used in the triangulation that it integrates on. - * - * @ingroup hp - */ - template - class FEValuesMap - { - public: - /** - * Make destructor virtual, - * since we have virtual - * functions in this class. - */ - virtual ~FEValuesMap (); - - /** - * Return a reference to the - * @p FEValues object selected - * by the last call to - * @p select_fe_values. The - * returned value is a constant - * reference since the only - * non-const function in the - * underlying class would be - * @p reinit, which you must - * not call directly any way; - * rather, use the @p reinit - * function of the - * hpFE*Values class. - */ - const FEValues & get_present_fe_values () const; - - protected: - /** - * Create an object of type - * @p FEValues for this - * particular finite - * element. Since the type of - * @p FEValues is not known - * here, we have to leave this - * to derived classes. - */ - virtual - FEValues * - create_fe_values (const FiniteElement &fe, - const unsigned int active_fe_index) const = 0; - - /** - * Select the @p FEValues - * object corresponding to the - * finite element object given - * as argument. If there is no - * such @p FEValues object - * yet, then one is created by - * calling @p create_fe_values - * with this finite element - * object. - * - * A non-const reference to - * this object is returned. - */ - FEValues & - select_fe_values (const FiniteElement &fe, - const unsigned int active_fe_index); - - - /** - * This field remembers the - * @p{active_fe_index} of the - * cell, which was used for - * the last call to @p{reinit}. - */ - unsigned int present_fe_index; - - private: - /** - * Have a map from pointers to - * finite elements to objects - * of type @p FEValues. - * - * Note that the shared pointer - * will make sure that the - * memory previously allocated - * is released when the - * destructor of this class is - * run. - */ - std::map >, unsigned int>, - boost::shared_ptr > fe_to_fe_values_map; - - /** - * Pointer to the @p FEValues - * object that is used for the - * present cell. This always - * points to one of the objects - * in the map above, and to - * which one is determined by - * the last call to - * @p select_fe_values. - */ - boost::shared_ptr present_fe_values; - - }; - - -/** - * Base class for the hp::FE*Values classes, storing the data that - * is common to them. The first template parameter denotes the space - * dimension we are in, the second the dimensionality of the object - * that we integrate on, i.e. for usual @p hp::FEValues it is equal to - * the first one, while for face integration it is one less. + * Base class for the hp::FE*Values classes, storing the data + * that is common to them. The first template parameter denotes the + * space dimension we are in, the second the dimensionality of the + * object that we integrate on, i.e. for usual @p hp::FEValues it is + * equal to the first one, while for face integration it is one + * less. The third template parameter indicates the type of underlying + * non-hp FE*Values base type, i.e. it could either be ::FEValues, + * ::FEFaceValues, or ::FESubfaceValues. * * @ingroup hp * * @author Wolfgang Bangerth, 2003 */ - template + template class FEValuesBase { - public: + public: /** * Constructor. Set the fields * of this class to the values @@ -163,10 +56,9 @@ namespace internal * to the constructor. */ FEValuesBase (const ::hp::MappingCollection &mapping_collection, + const ::hp::FECollection &fe_collection, const ::hp::QCollection &q_collection, const UpdateFlags update_flags); - - /** * Constructor. Set the fields * of this class to the values @@ -176,18 +68,59 @@ namespace internal * object for the mapping * object. */ - FEValuesBase (const ::hp::QCollection &q_collection, + FEValuesBase (const ::hp::FECollection &fe_collection, + const ::hp::QCollection &q_collection, const UpdateFlags update_flags); + + /** + * Return a reference to the + * @p FEValues object + * selected by the last call + * to + * select_fe_values(). select_fe_values() + * in turn is called when you + * called the @p reinit + * function of the + * hp::FE*Values + * class the last time. + */ + const FEValues & get_present_fe_values () const; + protected: + + /** + * Select a FEValues object + * suitable for the given FE, + * quadrature, and mapping + * indices. If such an object + * doesn't yet exist, create + * one. + * + * The function returns a + * writable reference so that + * derived classes can also + * reinit() the selected + * FEValues object. + */ + FEValues & + select_fe_values (const unsigned int fe_index, + const unsigned int mapping_index, + const unsigned int q_index); protected: /** - * A copy of the hp::MappingCollection - * object, which was specified - * upon construction of the - * object. + * A pointer to the + * collection of finite + * elements to be used. + */ + const SmartPointer > fe_collection; + + /** + * A pointer to the + * collection of mappings to + * be used. */ - SmartPointer > mapping_collection; + const SmartPointer > mapping_collection; /** * Copy of the quadrature @@ -197,6 +130,42 @@ namespace internal */ const ::hp::QCollection q_collection; + private: + /** + * A table in which we store + * pointers to fe_values + * objects for different + * finite element, mapping, + * and quadrature objects + * from our collection. The + * first index indicates the + * index of the finite + * element within the + * fe_collection, the second + * the index of the mapping + * within the mapping + * collection, and the last + * one the index of the + * quadrature formula within + * the q_collection. + * + * Initially, all entries + * have zero pointers, and we + * will allocate them lazily + * as needed in + * select_fe_values(). + */ + Table<3,boost::shared_ptr > fe_values_table; + + /** + * Set of indices pointing at + * the fe_values object + * selected last time the + * select_fe_value() function + * was called. + */ + TableIndices<3> present_fe_values_index; + /** * Values of the update flags as * given to the constructor. @@ -217,8 +186,7 @@ namespace hp * @ingroup hp */ template - class FEValues : public internal::hp::FEValuesMap >, - protected internal::hp::FEValuesBase + class FEValues : public internal::hp::FEValuesBase > { public: /** @@ -276,26 +244,108 @@ namespace hp /** * Reinitialize the object for - * the given cell. This selects - * the right @p FEValues object - * for the finite element in use - * by the cell given, and calling - * the @p reinit function on - * this object. + * the given cell. + * + * After the call, you can get + * an FEValues object using the + * get_present_fe_values() + * function that corresponds to + * the present cell. For this + * FEValues object, we use the + * additional arguments + * described below to determine + * which finite element, + * mapping, and quadrature + * formula to use. They are + * order in such a way that the + * arguments one may want to + * change most frequently come + * first. The rules for these + * arguments are as follows: + * + * If the @p fe_index argument + * to this function is left at + * its default value, then we + * use that finite element + * within the hp::FECollection + * passed to the constructor of + * this class with index given + * by + * cell-@>active_fe_index(). Consequently, + * the hp::FECollection + * argument given to this + * object should really be the + * same as that used in the + * construction of the + * hp::DofHandler associated + * with the present cell. On + * the other hand, if a value + * is given for this argument, + * it overrides the choice of + * cell-@>active_fe_index(). + * + * If the @p q_index argument + * is left at its default + * value, then we use that + * quadrature formula within + * the hp::QCollection passed + * to the constructor of this + * class with index given by + * cell-@>active_fe_index(), + * i.e. the same index as that + * of the finite element. In + * this case, there should be a + * corresponding quadrature + * formula for each finite + * element in the + * hp::FECollection. As a + * special case, if the + * quadrature collection + * contains only a single + * element (a frequent case if + * one wants to use the same + * quadrature object for all + * finite elements in an hp + * discretization, even if that + * may not be the most + * efficient), then this single + * quadrature is used unless a + * different value for this + * argument is specified. On + * the other hand, if a value + * is given for this argument, + * it overrides the choice of + * cell-@>active_fe_index() + * or the choice for the single + * quadrature. + * + * If the @p mapping_index + * argument is left at its + * default value, then we use + * that mapping object within + * the hp::MappingCollection + * passed to the constructor of + * this class with index given + * by + * cell-@>active_fe_index(), + * i.e. the same index as that + * of the finite + * element. As above, if the + * mapping collection contains + * only a single element (a + * frequent case if one wants + * to use a MappingQ1 object + * for all finite elements in + * an hp discretization), then + * this single mapping is used + * unless a different value for + * this argument is specified. */ void - reinit (const typename hp::DoFHandler::cell_iterator &cell); - - protected: - /** - * Create an object of type - * @p FEValues for this - * particular finite element. - */ - virtual - ::FEValues * - create_fe_values (const FiniteElement &fe, - const unsigned int active_fe_index) const; + reinit (const typename hp::DoFHandler::cell_iterator &cell, + const unsigned int q_index = deal_II_numbers::invalid_unsigned_int, + const unsigned int mapping_index = deal_II_numbers::invalid_unsigned_int, + const unsigned int fe_index = deal_II_numbers::invalid_unsigned_int); }; @@ -305,8 +355,7 @@ namespace hp * @ingroup hp */ template - class FEFaceValues : public internal::hp::FEValuesMap >, - protected internal::hp::FEValuesBase + class FEFaceValues : public internal::hp::FEValuesBase > { public: /** @@ -358,57 +407,114 @@ namespace hp * function. */ FEFaceValues (const hp::FECollection &fe_collection, - const hp::QCollection &q_collection, + const hp::QCollection &q_collection, const UpdateFlags update_flags); - - /** - * Reinitialize the object for - * the given cell. This selects - * the right @p FEValues object - * for the finite element in use - * by the cell given, and calling - * the @p reinit function on - * this object. - */ - void - reinit (const typename hp::DoFHandler::cell_iterator &cell, - const unsigned int face_no); - - /** * Reinitialize the object for - * the given cell. In this - * method, the user can specify - * which active_fe_index - * should be used to initialise - * the underlying FEValues - * object. This functionality - * is required, if the face terms - * between two cells with different - * polynomial degree should be - * assembled. In this case the - * values on the face of the - * lower order element have to - * be evaluated at the quadratrure - * points of the higher order - * element. + * the given cell and face. + * + * After the call, you can get + * an FEFaceValues object using the + * get_present_fe_values() + * function that corresponds to + * the present cell. For this + * FEFaceValues object, we use the + * additional arguments + * described below to determine + * which finite element, + * mapping, and quadrature + * formula to use. They are + * order in such a way that the + * arguments one may want to + * change most frequently come + * first. The rules for these + * arguments are as follows: + * + * If the @p fe_index argument + * to this function is left at + * its default value, then we + * use that finite element + * within the hp::FECollection + * passed to the constructor of + * this class with index given + * by + * cell-@>active_fe_index(). Consequently, + * the hp::FECollection + * argument given to this + * object should really be the + * same as that used in the + * construction of the + * hp::DofHandler associated + * with the present cell. On + * the other hand, if a value + * is given for this argument, + * it overrides the choice of + * cell-@>active_fe_index(). + * + * If the @p q_index argument + * is left at its default + * value, then we use that + * quadrature formula within + * the hp::QCollection passed + * to the constructor of this + * class with index given by + * cell-@>active_fe_index(), + * i.e. the same index as that + * of the finite element. In + * this case, there should be a + * corresponding quadrature + * formula for each finite + * element in the + * hp::FECollection. As a + * special case, if the + * quadrature collection + * contains only a single + * element (a frequent case if + * one wants to use the same + * quadrature object for all + * finite elements in an hp + * discretization, even if that + * may not be the most + * efficient), then this single + * quadrature is used unless a + * different value for this + * argument is specified. On + * the other hand, if a value + * is given for this argument, + * it overrides the choice of + * cell-@>active_fe_index() + * or the choice for the single + * quadrature. + * + * If the @p mapping_index + * argument is left at its + * default value, then we use + * that mapping object within + * the hp::MappingCollection + * passed to the constructor of + * this class with index given + * by + * cell-@>active_fe_index(), + * i.e. the same index as that + * of the finite + * element. As above, if the + * mapping collection contains + * only a single element (a + * frequent case if one wants + * to use a MappingQ1 object + * for all finite elements in + * an hp discretization), then + * this single mapping is used + * unless a different value for + * this argument is specified. */ void reinit (const typename hp::DoFHandler::cell_iterator &cell, const unsigned int face_no, - const unsigned int active_fe_index); - - protected: - /** - * Create an object of type - * @p FEValues for this - * particular finite element. - */ - virtual - ::FEFaceValues * - create_fe_values (const FiniteElement &fe, - const unsigned int active_fe_index) const; + const unsigned int q_index = deal_II_numbers::invalid_unsigned_int, + const unsigned int mapping_index = deal_II_numbers::invalid_unsigned_int, + const unsigned int fe_index = deal_II_numbers::invalid_unsigned_int); }; @@ -418,8 +524,7 @@ namespace hp * @ingroup hp */ template - class FESubfaceValues : public internal::hp::FEValuesMap >, - protected internal::hp::FEValuesBase + class FESubfaceValues : public internal::hp::FEValuesBase > { public: /** @@ -474,57 +579,91 @@ namespace hp const hp::QCollection &q_collection, const UpdateFlags update_flags); - /** * Reinitialize the object for - * the given cell. This selects - * the right @p FEValues object - * for the finite element in use - * by the cell given, and calling - * the @p reinit function on - * this object. - */ - void - reinit (const typename hp::DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no); - - - /** - * Reinitialize the object for - * the given cell. In this - * method, the user can specify - * which active_fe_index - * should be used to initialise - * the underlying FEValues - * object. This functionality - * is required, if the face terms - * between two cells with different - * polynomial degree should be - * assembled. In this case the - * values on the face of the - * lower order element have to - * be evaluated at the quadratrure - * points of the higher order - * element. + * the given cell, face, and subface. + * + * After the call, you can get + * an FESubfaceValues object using the + * get_present_fe_values() + * function that corresponds to + * the present cell. For this + * FESubfaceValues object, we use the + * additional arguments + * described below to determine + * which finite element, + * mapping, and quadrature + * formula to use. They are + * order in such a way that the + * arguments one may want to + * change most frequently come + * first. The rules for these + * arguments are as follows: + * + * If the @p q_index argument + * is left at its default + * value, then we use that + * quadrature formula within + * the hp::QCollection passed + * to the constructor of this + * class with index given by + * cell-@>active_fe_index(), + * i.e. the same index as that + * of the finite element. In + * this case, there should be a + * corresponding quadrature + * formula for each finite + * element in the + * hp::FECollection. As a + * special case, if the + * quadrature collection + * contains only a single + * element (a frequent case if + * one wants to use the same + * quadrature object for all + * finite elements in an hp + * discretization, even if that + * may not be the most + * efficient), then this single + * quadrature is used unless a + * different value for this + * argument is specified. On + * the other hand, if a value + * is given for this argument, + * it overrides the choice of + * cell-@>active_fe_index() + * or the choice for the single + * quadrature. + * + * If the @p mapping_index + * argument is left at its + * default value, then we use + * that mapping object within + * the hp::MappingCollection + * passed to the constructor of + * this class with index given + * by + * cell-@>active_fe_index(), + * i.e. the same index as that + * of the finite + * element. As above, if the + * mapping collection contains + * only a single element (a + * frequent case if one wants + * to use a MappingQ1 object + * for all finite elements in + * an hp discretization), then + * this single mapping is used + * unless a different value for + * this argument is specified. */ void reinit (const typename hp::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, - const unsigned int active_fe_index); - - - protected: - /** - * Create an object of type - * @p FEValues for this - * particular finite element. - */ - virtual - ::FESubfaceValues * - create_fe_values (const FiniteElement &fe, - const unsigned int active_fe_index) const; + const unsigned int q_index = deal_II_numbers::invalid_unsigned_int, + const unsigned int mapping_index = deal_II_numbers::invalid_unsigned_int, + const unsigned int fe_index = deal_II_numbers::invalid_unsigned_int); }; } @@ -536,12 +675,12 @@ namespace internal { namespace hp { - template + template inline const FEValues & - FEValuesMap::get_present_fe_values () const + FEValuesBase::get_present_fe_values () const { - return *present_fe_values; + return *fe_values_table(present_fe_values_index); } } diff --git a/deal.II/deal.II/source/fe/hp_fe_values.cc b/deal.II/deal.II/source/fe/hp_fe_values.cc index ad8fc8b421..9d38944f0e 100644 --- a/deal.II/deal.II/source/fe/hp_fe_values.cc +++ b/deal.II/deal.II/source/fe/hp_fe_values.cc @@ -16,65 +16,90 @@ #include -// -------------------------- FEValuesMap ------------------------- - namespace internal { namespace hp { - template - FEValuesMap::~FEValuesMap () - {} - - - template - FEValues & - FEValuesMap::select_fe_values (const FiniteElement &fe, - const unsigned int active_fe_index) - { - std::pair >, unsigned int> fe_pair= - std::make_pair(&fe, active_fe_index); - // check if the finite element - // does not exist as a key in the - // map - if (fe_to_fe_values_map.find (fe_pair) == fe_to_fe_values_map.end()) - // a-ha! doesn't yet, so let's - // make it up - fe_to_fe_values_map[fe_pair] - = boost::shared_ptr (create_fe_values (fe, active_fe_index)); - - - // now there definitely is one! - present_fe_values = fe_to_fe_values_map[fe_pair]; - - return *present_fe_values; - } - - // -------------------------- FEValuesBase ------------------------- - template - FEValuesBase::FEValuesBase ( - const ::hp::MappingCollection &mapping_collection, - const ::hp::QCollection &q_collection, - const UpdateFlags update_flags) + template + FEValuesBase:: + FEValuesBase (const ::hp::MappingCollection &mapping_collection, + const ::hp::FECollection &fe_collection, + const ::hp::QCollection &q_collection, + const UpdateFlags update_flags) : + fe_collection (&fe_collection), mapping_collection (&mapping_collection), q_collection (q_collection), + fe_values_table (fe_collection.size(), + mapping_collection.size(), + q_collection.size()), + present_fe_values_index (deal_II_numbers::invalid_unsigned_int, + deal_II_numbers::invalid_unsigned_int, + deal_II_numbers::invalid_unsigned_int), update_flags (update_flags) {} - template - FEValuesBase::FEValuesBase (const ::hp::QCollection &q_collection, - const UpdateFlags update_flags) + template + FEValuesBase:: + FEValuesBase (const ::hp::FECollection &fe_collection, + const ::hp::QCollection &q_collection, + const UpdateFlags update_flags) : + fe_collection (&fe_collection), mapping_collection (&::hp::StaticMappingQ1::mapping_collection), q_collection (q_collection), + fe_values_table (fe_collection.size(), + 1, + q_collection.size()), + present_fe_values_index (deal_II_numbers::invalid_unsigned_int, + deal_II_numbers::invalid_unsigned_int, + deal_II_numbers::invalid_unsigned_int), update_flags (update_flags) {} - + + + + template + FEValues & + FEValuesBase:: + select_fe_values (const unsigned int fe_index, + const unsigned int mapping_index, + const unsigned int q_index) + { + Assert (fe_index < fe_collection->size(), + ExcIndexRange (fe_index, 0, fe_collection->size())); + Assert (mapping_index < mapping_collection->size(), + ExcIndexRange (mapping_index, 0, mapping_collection->size())); + Assert (q_index < q_collection.size(), + ExcIndexRange (q_index, 0, q_collection.size())); + + + // set the triple of indices + // that we want to work with + present_fe_values_index = TableIndices<3> (fe_index, + mapping_index, + q_index); + + // first check whether we + // already have an object for + // this particular combination + // of indices + if (fe_values_table(present_fe_values_index) == 0) + fe_values_table(present_fe_values_index) + = + boost::shared_ptr + (new FEValues ((*mapping_collection)[mapping_index], + (*fe_collection)[fe_index], + q_collection[q_index], + update_flags)); + + // now there definitely is one! + return *fe_values_table(present_fe_values_index); + } } } @@ -88,44 +113,70 @@ namespace hp template FEValues::FEValues (const hp::MappingCollection &mapping, - const hp::FECollection &/*fe_collection*/, + const hp::FECollection &fe_collection, const hp::QCollection &q_collection, - const UpdateFlags update_flags) + const UpdateFlags update_flags) : - internal::hp::FEValuesBase (mapping, - q_collection, - update_flags) + internal::hp::FEValuesBase > (mapping, + fe_collection, + q_collection, + update_flags) {} template - FEValues::FEValues (const hp::FECollection &/*fe_collection*/, + FEValues::FEValues (const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags) : - internal::hp::FEValuesBase (q_collection, + internal::hp::FEValuesBase > (fe_collection, + q_collection, update_flags) {} template void - FEValues::reinit (const typename hp::DoFHandler::cell_iterator &cell) - { - this->present_fe_index = cell->active_fe_index (); - this->select_fe_values (cell->get_fe(), this->present_fe_index).reinit (cell); - } - - - - template - ::FEValues * - FEValues::create_fe_values (const FiniteElement &fe, - const unsigned int active_fe_index) const + FEValues::reinit (const typename hp::DoFHandler::cell_iterator &cell, + const unsigned int q_index, + const unsigned int mapping_index, + const unsigned int fe_index) { - return new ::FEValues ( - (*this->mapping_collection)[active_fe_index], fe, - this->q_collection[active_fe_index], this->update_flags); + // determine which indices we + // should actually use + unsigned int real_q_index = q_index, + real_mapping_index = mapping_index, + real_fe_index = fe_index; + + if (real_q_index == deal_II_numbers::invalid_unsigned_int) + if (q_collection.size() > 1) + real_q_index = cell->active_fe_index(); + else + real_q_index = 0; + + if (real_mapping_index == deal_II_numbers::invalid_unsigned_int) + if (mapping_collection->size() > 1) + real_mapping_index = cell->active_fe_index(); + else + real_mapping_index = 0; + + if (real_fe_index == deal_II_numbers::invalid_unsigned_int) + real_fe_index = cell->active_fe_index(); + + // some checks + Assert (real_q_index < q_collection.size(), + ExcIndexRange (real_q_index, 0, q_collection.size())); + Assert (real_mapping_index < mapping_collection->size(), + ExcIndexRange (real_mapping_index, 0, mapping_collection->size())); + Assert (real_fe_index < fe_collection->size(), + ExcIndexRange (real_fe_index, 0, fe_collection->size())); + + // now finally actually get the + // corresponding object and + // initialize it + this->select_fe_values (real_fe_index, + real_mapping_index, + real_q_index).reinit (cell); } @@ -134,55 +185,71 @@ namespace hp template FEFaceValues::FEFaceValues (const hp::MappingCollection &mapping, - const hp::FECollection &/*fe_collection*/, + const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags) : - internal::hp::FEValuesBase (mapping, + internal::hp::FEValuesBase > (mapping, + fe_collection, q_collection, update_flags) {} template - FEFaceValues::FEFaceValues (const hp::FECollection &/*fe_collection*/, + FEFaceValues::FEFaceValues (const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags) : - internal::hp::FEValuesBase (q_collection, + internal::hp::FEValuesBase > (fe_collection, + q_collection, update_flags) {} - template - void - FEFaceValues::reinit (const typename hp::DoFHandler::cell_iterator &cell, - const unsigned int face_no) - { - this->present_fe_index = cell->active_fe_index (); - this->select_fe_values (cell->get_fe(), this->present_fe_index).reinit (cell, face_no); - } - - template void FEFaceValues::reinit (const typename hp::DoFHandler::cell_iterator &cell, const unsigned int face_no, - const unsigned int active_fe_index) - { - this->present_fe_index = active_fe_index; - this->select_fe_values (cell->get_fe(), active_fe_index).reinit (cell, face_no); - } - - - template - ::FEFaceValues * - FEFaceValues::create_fe_values (const FiniteElement &fe, - const unsigned int active_fe_index) const + const unsigned int q_index, + const unsigned int mapping_index, + const unsigned int fe_index) { - return new ::FEFaceValues ( - (*this->mapping_collection)[active_fe_index], fe, - this->q_collection[active_fe_index], this->update_flags); + // determine which indices we + // should actually use + unsigned int real_q_index = q_index, + real_mapping_index = mapping_index, + real_fe_index = fe_index; + + if (real_q_index == deal_II_numbers::invalid_unsigned_int) + if (q_collection.size() > 1) + real_q_index = cell->active_fe_index(); + else + real_q_index = 0; + + if (real_mapping_index == deal_II_numbers::invalid_unsigned_int) + if (mapping_collection->size() > 1) + real_mapping_index = cell->active_fe_index(); + else + real_mapping_index = 0; + + if (real_fe_index == deal_II_numbers::invalid_unsigned_int) + real_fe_index = cell->active_fe_index(); + + // some checks + Assert (real_q_index < q_collection.size(), + ExcIndexRange (real_q_index, 0, q_collection.size())); + Assert (real_mapping_index < mapping_collection->size(), + ExcIndexRange (real_mapping_index, 0, mapping_collection->size())); + Assert (real_fe_index < fe_collection->size(), + ExcIndexRange (real_fe_index, 0, fe_collection->size())); + + // now finally actually get the + // corresponding object and + // initialize it + this->select_fe_values (real_fe_index, + real_mapping_index, + real_q_index).reinit (cell, face_no); } @@ -191,57 +258,72 @@ namespace hp template FESubfaceValues::FESubfaceValues (const hp::MappingCollection &mapping, - const hp::FECollection &/*fe_collection*/, + const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags) : - internal::hp::FEValuesBase (mapping, + internal::hp::FEValuesBase > (mapping, + fe_collection, q_collection, update_flags) {} template - FESubfaceValues::FESubfaceValues (const hp::FECollection &/*fe_collection*/, + FESubfaceValues::FESubfaceValues (const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags) : - internal::hp::FEValuesBase (q_collection, + internal::hp::FEValuesBase > (fe_collection, + q_collection, update_flags) {} - template - void - FESubfaceValues::reinit (const typename hp::DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no) - { - this->present_fe_index = cell->active_fe_index (); - this->select_fe_values (cell->get_fe(), this->present_fe_index).reinit (cell, face_no, subface_no); - } - - template void FESubfaceValues::reinit (const typename hp::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, - const unsigned int active_fe_index) + const unsigned int q_index, + const unsigned int mapping_index, + const unsigned int fe_index) { - this->present_fe_index = active_fe_index; - this->select_fe_values (cell->get_fe(), active_fe_index).reinit (cell, face_no, subface_no); - } - - - template - ::FESubfaceValues * - FESubfaceValues::create_fe_values (const FiniteElement &fe, - const unsigned int active_fe_index) const - { - return new ::FESubfaceValues ( - (*this->mapping_collection)[active_fe_index], fe, - this->q_collection[active_fe_index], this->update_flags); + // determine which indices we + // should actually use + unsigned int real_q_index = q_index, + real_mapping_index = mapping_index, + real_fe_index = fe_index; + + if (real_q_index == deal_II_numbers::invalid_unsigned_int) + if (q_collection.size() > 1) + real_q_index = cell->active_fe_index(); + else + real_q_index = 0; + + if (real_mapping_index == deal_II_numbers::invalid_unsigned_int) + if (mapping_collection->size() > 1) + real_mapping_index = cell->active_fe_index(); + else + real_mapping_index = 0; + + if (real_fe_index == deal_II_numbers::invalid_unsigned_int) + real_fe_index = cell->active_fe_index(); + + // some checks + Assert (real_q_index < q_collection.size(), + ExcIndexRange (real_q_index, 0, q_collection.size())); + Assert (real_mapping_index < mapping_collection->size(), + ExcIndexRange (real_mapping_index, 0, mapping_collection->size())); + Assert (real_fe_index < fe_collection->size(), + ExcIndexRange (real_fe_index, 0, fe_collection->size())); + + // now finally actually get the + // corresponding object and + // initialize it + this->select_fe_values (real_fe_index, + real_mapping_index, + real_q_index).reinit (cell, face_no, subface_no); } } @@ -251,9 +333,13 @@ namespace internal { namespace hp { - template class FEValuesBase; + template class FEValuesBase >; #if deal_II_dimension >= 2 - template class FEValuesBase; + template class FEValuesBase >; + template class FEValuesBase >; #endif } } @@ -266,10 +352,3 @@ namespace hp template class FESubfaceValues; #endif } - -// Putting the following explicit instantiations into the brackets -// of the appropriate namespace somehow causes problems with the -// Apple gcc3.3. Therefore these are separated. -template class internal::hp::FEValuesMap >; -template class internal::hp::FEValuesMap >; -template class internal::hp::FEValuesMap >; diff --git a/deal.II/examples/step-21/step-21.cc b/deal.II/examples/step-21/step-21.cc index 564f3d787b..0e4a21fd0d 100644 --- a/deal.II/examples/step-21/step-21.cc +++ b/deal.II/examples/step-21/step-21.cc @@ -1010,9 +1010,9 @@ void DGMethod::assemble_system1 () // the face quadrature rule of the // higher order element // will be used. - const unsigned int use_fe_index = - neighbor_child->active_fe_index () > cell->active_fe_index () ? - neighbor_child->active_fe_index () : cell->active_fe_index (); + const unsigned int quadrature_index = + std::max (neighbor_child->active_fe_index (), + cell->active_fe_index ()); // As these are @@ -1067,8 +1067,8 @@ void DGMethod::assemble_system1 () // of the // neighboring // child cell. - fe_v_subface_x.reinit (cell, face_no, subface_no, use_fe_index); - fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, use_fe_index); + fe_v_subface_x.reinit (cell, face_no, subface_no, quadrature_index); + fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, quadrature_index); dg.assemble_face_term1(fe_v_subface_x.get_present_fe_values (), fe_v_face_neighbor_x.get_present_fe_values (), @@ -1124,9 +1124,9 @@ void DGMethod::assemble_system1 () // quadrature rule // of higher order // cell. - const unsigned int use_fe_index = - neighbor->active_fe_index () > cell->active_fe_index () ? - neighbor->active_fe_index () : cell->active_fe_index (); + const unsigned int quadrature_index = + std::max (neighbor->active_fe_index (), + cell->active_fe_index ()); // We reinit // the @@ -1140,8 +1140,8 @@ void DGMethod::assemble_system1 () // the // corresponding // face terms. - fe_v_face_x.reinit (cell, face_no, use_fe_index); - fe_v_face_neighbor_x.reinit (neighbor, neighbor2, use_fe_index); + fe_v_face_x.reinit (cell, face_no, quadrature_index); + fe_v_face_neighbor_x.reinit (neighbor, neighbor2, quadrature_index); dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (), fe_v_face_neighbor_x.get_present_fe_values (), @@ -1193,9 +1193,9 @@ void DGMethod::assemble_system1 () // quadrature rule // of higher order // cell. - const unsigned int use_fe_index = - neighbor->active_fe_index () > cell->active_fe_index () ? - neighbor->active_fe_index () : cell->active_fe_index (); + const unsigned int quadrature_index = + std::max (neighbor->active_fe_index (), + cell->active_fe_index ()); // Reinit the // appropriate @@ -1203,9 +1203,9 @@ void DGMethod::assemble_system1 () // and assemble // the face // terms. - fe_v_face_x.reinit (cell, face_no, use_fe_index); + fe_v_face_x.reinit (cell, face_no, quadrature_index); fe_v_subface_neighbor_x.reinit (neighbor, neighbor_face_no, - neighbor_subface_no, use_fe_index); + neighbor_subface_no, quadrature_index); dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (), fe_v_subface_neighbor_x.get_present_fe_values (), @@ -1390,9 +1390,9 @@ void DGMethod::assemble_system2 () neighbor2,subface_no)); const unsigned int dofs_on_neighbor_child = neighbor_child->get_fe().dofs_per_cell; - const unsigned int use_fe_index = - neighbor_child->active_fe_index () > cell->active_fe_index () ? - neighbor_child->active_fe_index () : cell->active_fe_index (); + const unsigned int quadrature_index = + std::max (neighbor_child->active_fe_index (), + cell->active_fe_index ()); Assert (neighbor_child->face(neighbor2) == face->child(subface_no), ExcInternalError()); @@ -1402,8 +1402,8 @@ void DGMethod::assemble_system2 () u_vn_matrix = 0; un_vn_matrix = 0; - fe_v_subface_x.reinit (cell, face_no, subface_no, use_fe_index); - fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, use_fe_index); + fe_v_subface_x.reinit (cell, face_no, subface_no, quadrature_index); + fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, quadrature_index); dg.assemble_face_term2(fe_v_subface_x.get_present_fe_values (), fe_v_face_neighbor_x.get_present_fe_values (), @@ -1441,16 +1441,16 @@ void DGMethod::assemble_system2 () { const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); - const unsigned int use_fe_index = - neighbor->active_fe_index () > cell->active_fe_index () ? - neighbor->active_fe_index () : cell->active_fe_index (); + const unsigned int quadrature_index = + std::max (neighbor->active_fe_index (), + cell->active_fe_index ()); un_v_matrix = 0; u_vn_matrix = 0; un_vn_matrix = 0; - fe_v_face_x.reinit (cell, face_no, use_fe_index); - fe_v_face_neighbor_x.reinit (neighbor, neighbor2, use_fe_index); + fe_v_face_x.reinit (cell, face_no, quadrature_index); + fe_v_face_neighbor_x.reinit (neighbor, neighbor2, quadrature_index); dg.assemble_face_term2(fe_v_face_x.get_present_fe_values (), fe_v_face_neighbor_x.get_present_fe_values (), @@ -1658,7 +1658,7 @@ void DGMethod::output_results (const unsigned int cycle) const Assert (cycle < 10, ExcInternalError()); // filename += ".gnuplot"; - filename += ".vtk"; + filename += ".gmv"; deallog << "Writing solution to <" << filename << ">..." << std::endl << std::endl; std::ofstream gnuplot_output (filename.c_str()); @@ -1670,7 +1670,7 @@ void DGMethod::output_results (const unsigned int cycle) const data_out.build_patches (5); // data_out.write_gnuplot(gnuplot_output); - data_out.write_vtk(gnuplot_output); + data_out.write_gmv(gnuplot_output); }