From: wolf Date: Sun, 27 Jun 2004 21:23:35 +0000 (+0000) Subject: Merge most of the changes in the fe subdirectory from the hp branch to mainline,... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c7117dcc3737dd460db820967ade312edaff2269;p=dealii-svn.git Merge most of the changes in the fe subdirectory from the hp branch to mainline, so as to avoid further divergence of the branches. git-svn-id: https://svn.dealii.org/trunk@9470 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 19779ee73a..e3c9e4f3e9 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -23,6 +23,7 @@ template class FEValues; template class FEFaceValues; template class FESubfaceValues; template class FESystem; +template class FECollection; /*!@addtogroup febase */ /*@{*/ @@ -242,7 +243,8 @@ class FiniteElement : public FiniteElementBase * a copy constructor. * * This function is needed by the - * constructors of @p FESystem. + * constructors of FESystem as well + * as by the FECollection class. */ virtual FiniteElement *clone() const = 0; @@ -309,7 +311,7 @@ class FiniteElement : public FiniteElementBase */ virtual void fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_internal, typename Mapping::InternalDataBase &fe_internal, @@ -328,7 +330,7 @@ class FiniteElement : public FiniteElementBase */ virtual void fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_internal, @@ -348,7 +350,7 @@ class FiniteElement : public FiniteElementBase */ virtual void fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, @@ -365,6 +367,7 @@ class FiniteElement : public FiniteElementBase friend class FEFaceValues; friend class FESubfaceValues; friend class FESystem; + friend class FECollection; }; /*@}*/ diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index 551f7058df..df31b9fddf 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -31,7 +31,6 @@ template class FESystem; /*!@addtogroup febase */ /*@{*/ - /** * Dimension independent data for finite elements. See the derived * class FiniteElementBase class for information on its use. All @@ -904,7 +903,7 @@ class FiniteElementBase : public Subscriptor, * @p system_to_base_index * function. */ - std::pair + std::pair system_to_component_index (const unsigned int index) const; /** @@ -912,7 +911,7 @@ class FiniteElementBase : public Subscriptor, * shape functions and their * indices on a face. */ - std::pair + std::pair face_system_to_component_index (const unsigned int index) const; /** @@ -1661,7 +1660,7 @@ class FiniteElementBase : public Subscriptor, * gradients. */ void compute_2nd (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int offset, typename Mapping::InternalDataBase &mapping_internal, InternalDataBase &fe_internal, diff --git a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h index 928f1c2426..191c368da1 100644 --- a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h +++ b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h @@ -310,7 +310,7 @@ class FE_DGPNonparametric : public FiniteElement */ virtual void fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_internal, typename Mapping::InternalDataBase &fe_internal, @@ -323,7 +323,7 @@ class FE_DGPNonparametric : public FiniteElement */ virtual void fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_internal, @@ -337,7 +337,7 @@ class FE_DGPNonparametric : public FiniteElement */ virtual void fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, diff --git a/deal.II/deal.II/include/fe/fe_nedelec.h b/deal.II/deal.II/include/fe/fe_nedelec.h index 4755163d5a..330da4aae8 100644 --- a/deal.II/deal.II/include/fe/fe_nedelec.h +++ b/deal.II/deal.II/include/fe/fe_nedelec.h @@ -413,7 +413,7 @@ class FE_Nedelec : public FiniteElement */ virtual void fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_internal, typename Mapping::InternalDataBase &fe_internal, @@ -426,7 +426,7 @@ class FE_Nedelec : public FiniteElement */ virtual void fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_internal, @@ -440,7 +440,7 @@ class FE_Nedelec : public FiniteElement */ virtual void fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, diff --git a/deal.II/deal.II/include/fe/fe_poly.h b/deal.II/deal.II/include/fe/fe_poly.h index 39818dcba2..cc54a38438 100644 --- a/deal.II/deal.II/include/fe/fe_poly.h +++ b/deal.II/deal.II/include/fe/fe_poly.h @@ -225,7 +225,7 @@ class FE_Poly : public FiniteElement */ virtual void fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_internal, typename Mapping::InternalDataBase &fe_internal, @@ -237,7 +237,7 @@ class FE_Poly : public FiniteElement */ virtual void fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_internal, @@ -250,7 +250,7 @@ class FE_Poly : public FiniteElement */ virtual void fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, diff --git a/deal.II/deal.II/include/fe/fe_raviart_thomas.h b/deal.II/deal.II/include/fe/fe_raviart_thomas.h index 9532dc1b4d..e4d3c4a2c2 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -318,7 +318,7 @@ class FE_RaviartThomas : public FiniteElement */ virtual void fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_internal, typename Mapping::InternalDataBase &fe_internal, @@ -331,7 +331,7 @@ class FE_RaviartThomas : public FiniteElement */ virtual void fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_internal, @@ -345,7 +345,7 @@ class FE_RaviartThomas : public FiniteElement */ virtual void fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index a4e4b870e9..63d64dc911 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -448,7 +448,7 @@ class FESystem : public FiniteElement */ virtual void fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fe_data, @@ -466,7 +466,7 @@ class FESystem : public FiniteElement */ virtual void fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, @@ -485,7 +485,7 @@ class FESystem : public FiniteElement */ virtual void fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, @@ -518,7 +518,7 @@ class FESystem : public FiniteElement */ template void compute_fill (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index f08a705684..8072d3dbad 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -15,10 +15,16 @@ +#include +#include + +#include + template class FullMatrix; +template class Vector; template class FiniteElement; template class DoFHandler; -template class Vector; +template class hpDoFHandler; template class FE_Q; class ConstraintMatrix; @@ -31,6 +37,7 @@ class ConstraintMatrix; /*!@addtogroup febase */ /*@{*/ + /** * This class performs interpolations and extrapolations of discrete * functions of one @p FiniteElement @p fe1 to another @p FiniteElement @@ -41,7 +48,7 @@ class ConstraintMatrix; * $id-I_h$ that is needed for evaluating $(id-I_h)z$ for e.g. the * dual solution $z$. * - * @author Ralf Hartmann, Guido Kanschat, 2000, 2004 + * @author Ralf Hartmann, 2000; Wolfgang Bangerth, 2003, Guido Kanschat, 2000, 2004 */ class FETools { @@ -88,9 +95,11 @@ class FETools * is an embedding matrix. */ template - static void get_interpolation_matrix(const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &interpolation_matrix); + static + void + get_interpolation_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix); /** * Gives the interpolation matrix @@ -111,9 +120,11 @@ class FETools * only the unit matrix. */ template - static void get_back_interpolation_matrix(const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &interpolation_matrix); + static + void + get_back_interpolation_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix); /** * Gives the unit matrix minus the @@ -136,8 +147,8 @@ class FETools static void get_interpolation_difference_matrix(const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &difference_matrix); + const FiniteElement &fe2, + FullMatrix &difference_matrix); /** * Compute the local @@ -207,11 +218,16 @@ class FETools * your hanging node constraints * object. */ - template - static void interpolate (const DoFHandler& dof1, - const InVector& u1, - const DoFHandler& dof2, - OutVector& u2); + template class DH1, + template class DH2, + class InVector, class OutVector> + static + void + interpolate (const DH1 &dof1, + const InVector &u1, + const DH2 &dof2, + OutVector &u2); /** * Gives the interpolation of a @@ -252,6 +268,23 @@ class FETools const DoFHandler& dof2, const ConstraintMatrix& constraints, OutVector& u2); + + + /** + * Same as last function, except + * that one or both of the dof + * handler objects might be of + * type @p hpDoFHandler. + */ + template class DH1, + template class DH2, + class InVector, class OutVector> + static void interpolate (const DH1 &dof1, + const InVector &u1, + const DH2 &dof2, + const ConstraintMatrix &constraints, + OutVector& u2); /** @@ -281,12 +314,26 @@ class FETools * mapping. */ template - static void back_interpolate (const DoFHandler& dof1, - const InVector& u1, - const FiniteElement& fe2, - OutVector& u1_interpolated); + static void back_interpolate (const DoFHandler &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_interpolated); - /** + /** + * Same as last function, except + * that the dof handler objects + * might be of type + * @p hpDoFHandler. + */ + template class DH, + class InVector, class OutVector> + static void back_interpolate (const DH &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_interpolated); + + /** * Gives the interpolation of the * @p dof1-function @p u1 to a * @p dof2-function, and diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 0b1968b28c..18a620e161 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -21,6 +21,10 @@ #include #include #include +#include +#include +#include +#include #include #include #include @@ -28,11 +32,15 @@ #include #include #include +#include +#include #include +#include template class Quadrature; + /*!@addtogroup febase */ /*@{*/ @@ -307,7 +315,7 @@ class FEValuesData * updating. Still, it is recommended to give ALL needed update flags * to @p FEValues. * - * @author Wolfgang Bangerth, 1998, Guido Kanschat, 2001 + * @author Wolfgang Bangerth, 1998, 2003, Guido Kanschat, 2001 */ template class FEValuesBase : protected FEValuesData @@ -824,10 +832,10 @@ class FEValuesBase : protected FEValuesData const FiniteElement & get_fe () const; /** - * Constant reference to the - * current cell + * Return a triangulation + * iterator to the current cell. */ - typename Triangulation::cell_iterator get_cell () const; + const typename Triangulation::cell_iterator get_cell () const; /** * Determine an estimate for the @@ -880,6 +888,469 @@ class FEValuesBase : protected FEValuesData DeclException0 (ExcFENotPrimitive); protected: + /** + * Objects of the @p FEValues + * class need to store a pointer + * (i.e. an iterator) to the + * present cell in order to be + * able to extract the values of + * the degrees of freedom on this + * cell in the + * @p get_function_values and + * assorted functions. On the + * other hand, this class should + * also work for different + * iterators, as long as they + * have the same interface to + * extract the DoF values (i.e., + * for example, they need to have + * a @p get_interpolated_dof_values + * function). + * + * This calls for a common base + * class of iterator classes, and + * making the functions we need + * here @p virtual. On the other + * hand, this is the only place + * in the library where we need + * this, and introducing a base + * class of iterators and making + * a function virtual penalizes + * @em{all} users of the + * iterators, which are basically + * intended as very fast accessor + * functions. So we do not want + * to do this. Rather, what we do + * here is making the functions + * we need virtual only for use + * with @em{this class}. The idea + * is the following: have a + * common base class which + * declares some pure virtual + * functions, and for each + * possible iterator type, we + * have a derived class which + * stores the iterator to the + * cell and implements these + * functions. Since the iterator + * classes have the same + * interface, we can make the + * derived classes a template, + * templatized on the iterator + * type. + * + * This way, the use of virtual + * functions is restricted to + * only this class, and other + * users of iterators do not have + * to bear the negative effects. + * + * @author Wolfgang Bangerth, 2003 + */ + class CellIteratorBase + { + public: + /** + * Destructor. Made virtual + * since we store only + * pointers to the base + * class. + */ + virtual ~CellIteratorBase (); + + /** + * Conversion operator to an + * iterator for + * triangulations. This + * conversion is implicit for + * the original iterators, + * since they are derived + * classes. However, since + * here we have kind of a + * parallel class hierarchy, + * we have to have a + * conversion operator. + */ + virtual + operator const typename Triangulation::cell_iterator () const = 0; + + /** + * Return the number of + * degrees of freedom the DoF + * handler object has to + * which the iterator belongs + * to. + */ + virtual + unsigned int + n_dofs_for_dof_handler () const = 0; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const Vector &in, + Vector &out) const = 0; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const Vector &in, + Vector &out) const = 0; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const BlockVector &in, + Vector &out) const = 0; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const BlockVector &in, + Vector &out) const = 0; + +#ifdef DEAL_II_USE_PETSC + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const PETScWrappers::Vector &in, + Vector &out) const = 0; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const PETScWrappers::BlockVector &in, + Vector &out) const = 0; +#endif + }; + + + /** + * Implementation of derived + * classes of the + * CellIteratorBase + * interface. See there for a + * description of the use of + * these classes. + * + * @author Wolfgang Bangerth, 2003 + */ + template + class CellIterator : public CellIteratorBase + { + public: + /** + * Constructor. Take an + * iterator and store it in + * this class. + */ + CellIterator (const CI &cell); + + /** + * Conversion operator to an + * iterator for + * triangulations. This + * conversion is implicit for + * the original iterators, + * since they are derived + * classes. However, since + * here we have kind of a + * parallel class hierarchy, + * we have to have a + * conversion operator. + */ + virtual + operator const typename Triangulation::cell_iterator () const; + + /** + * Return the number of + * degrees of freedom the DoF + * handler object has to + * which the iterator belongs + * to. + */ + virtual + unsigned int + n_dofs_for_dof_handler () const; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const Vector &in, + Vector &out) const; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const Vector &in, + Vector &out) const; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const BlockVector &in, + Vector &out) const; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const BlockVector &in, + Vector &out) const; + +#ifdef DEAL_II_USE_PETSC + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const PETScWrappers::Vector &in, + Vector &out) const; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const PETScWrappers::BlockVector &in, + Vector &out) const; +#endif + private: + /** + * Copy of the iterator which + * we use in this object. + */ + const CI cell; + }; + + + /** + * Implementation of a derived + * class of the + * CellIteratorBase + * interface. See there for a + * description of the use of + * these classes. + * + * This class is basically a + * specialization of the general + * template for iterators into + * @p Triangulation objects (but + * since C++ does not allow + * something like this for nested + * classes, it runs under a + * separate name). Since these do + * not implement the interface + * that we would like to call, + * the functions of this class + * cannot be implemented + * meaningfully. However, most + * functions of the @p FEValues + * class do not make any use of + * degrees of freedom at all, so + * it should be possible to call + * @p FEValues::reinit with a tria + * iterator only; this class + * makes this possible, but + * whenever one of the functions + * of @p FEValues tries to call + * any of the functions of this + * class, an exception will be + * raised reminding the user that + * if she wants to use these + * features, then the + * @p FEValues object has to be + * reinitialized with a cell + * iterator that allows to + * extract degree of freedom + * information. + * + * @author Wolfgang Bangerth, 2003 + */ + class TriaCellIterator : public CellIteratorBase + { + public: + /** + * Constructor. Take an + * iterator and store it in + * this class. + */ + TriaCellIterator (const typename Triangulation::cell_iterator &cell); + + /** + * Conversion operator to an + * iterator for + * triangulations. This + * conversion is implicit for + * the original iterators, + * since they are derived + * classes. However, since + * here we have kind of a + * parallel class hierarchy, + * we have to have a + * conversion operator. Here, + * the conversion is trivial, + * from and to the same time. + */ + virtual + operator const typename Triangulation::cell_iterator () const; + + /** + * Implement the respective + * function of the base + * class. Since this is not + * possible, we just raise an + * error. + */ + virtual + unsigned int + n_dofs_for_dof_handler () const; + + /** + * Implement the respective + * function of the base + * class. Since this is not + * possible, we just raise an + * error. + */ + virtual + void + get_interpolated_dof_values (const Vector &in, + Vector &out) const; + + /** + * Implement the respective + * function of the base + * class. Since this is not + * possible, we just raise an + * error. + */ + virtual + void + get_interpolated_dof_values (const Vector &in, + Vector &out) const; + + /** + * Implement the respective + * function of the base + * class. Since this is not + * possible, we just raise an + * error. + */ + virtual + void + get_interpolated_dof_values (const BlockVector &in, + Vector &out) const; + + /** + * Implement the respective + * function of the base + * class. Since this is not + * possible, we just raise an + * error. + */ + virtual + void + get_interpolated_dof_values (const BlockVector &in, + Vector &out) const; + +#ifdef DEAL_II_USE_PETSC + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const PETScWrappers::Vector &in, + Vector &out) const; + + /** + * Call + * @p get_interpolated_dof_values + * of the iterator with the + * given arguments. + */ + virtual + void + get_interpolated_dof_values (const PETScWrappers::BlockVector &in, + Vector &out) const; +#endif + private: + /** + * Copy of the iterator which + * we use in this object. + */ + const typename Triangulation::cell_iterator cell; + + /** + * String to be displayed + * whenever one of the + * functions of this class is + * called. Make it a static + * member variable, since we + * show the same message for + * all member functions. + */ + static const char * const message_string; + }; + /** * Store the cell selected last time * the @p reinit function was called @@ -887,8 +1358,8 @@ class FEValuesBase : protected FEValuesData * to the get_function_* functions * safer. */ - typename DoFHandler::cell_iterator present_cell; - + std::auto_ptr present_cell; + /** * Storage for the mapping object. */ @@ -937,6 +1408,23 @@ class FEValuesBase : protected FEValuesData * @p MappingQ1 implicitly. */ static const Mapping &get_default_mapping(); + + private: + /** + * Copy constructor. Since + * objects of this class are not + * copyable, we make it private, + * and also do not implement it. + */ + FEValuesBase (const FEValuesBase &); + + /** + * Copy operator. Since + * objects of this class are not + * copyable, we make it private, + * and also do not implement it. + */ + FEValuesBase & operator= (const FEValuesBase &); }; @@ -975,11 +1463,62 @@ class FEValues : public FEValuesBase const UpdateFlags update_flags); /** - * Reinitialize the gradients, Jacobi - * determinants, etc for the given cell - * and the given finite element. + * Reinitialize the gradients, + * Jacobi determinants, etc for + * the given cell of type + * "iterator into a DoFHandler + * object", and the finite + * element associated with this + * object. It is assumed that the + * finite element used by the + * given cell is also the one + * used by this @p FEValues + * object. */ - void reinit (const typename DoFHandler::cell_iterator &); + void reinit (const typename DoFHandler::cell_iterator &cell); + + /** + * Reinitialize the gradients, + * Jacobi determinants, etc for + * the given cell of type + * "iterator into a MGDoFHandler + * object", and the finite + * element associated with this + * object. It is assumed that the + * finite element used by the + * given cell is also the one + * used by this @p FEValues + * object. + */ + void reinit (const typename MGDoFHandler::cell_iterator &cell); + + /** + * Reinitialize the gradients, + * Jacobi determinants, etc for + * the given cell of type + * "iterator into a Triangulation + * object", and the given finite + * element. Since iterators into + * triangulation alone only + * convey information about the + * geometry of a cell, but not + * about degrees of freedom + * possibly associated with this + * cell, you will not be able to + * call some functions of this + * class if they need information + * about degrees of + * freedom. These functions are, + * above all, the + * get_function_value/grads/2nd_derivatives + * functions. If you want to call + * these functions, you have to + * call the @p reinit variants + * that take iterators into + * @p DoFHandler or other DoF + * handler type objects. + */ + void reinit (const typename Triangulation::cell_iterator &cell); /** * Return a reference to the copy @@ -995,6 +1534,13 @@ class FEValues : public FEValuesBase */ unsigned int memory_consumption () const; + /** + * Return a reference to this + * very object. + */ +//TODO: explain reason + const FEValues & get_present_fe_values () const; + private: /** * Store a copy of the quadrature @@ -1007,6 +1553,20 @@ class FEValues : public FEValuesBase * constructors. */ void initialize (const UpdateFlags update_flags); + + /** + * The reinit() functions do + * only that part of the work + * that requires knowledge of the + * type of iterator. After + * setting present_cell(), + * they pass on to this function, + * which does the real work, and + * which is independent of the + * actual type of the cell + * iterator. + */ + void do_reinit (); }; @@ -1139,18 +1699,18 @@ class FEFaceValues : public FEFaceValuesBase * matching the quadrature rule * and update flags. */ - FEFaceValues (const Mapping& mapping, - const FiniteElement& fe, - const Quadrature& quadrature, - const UpdateFlags update_flags); + FEFaceValues (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature &quadrature, + const UpdateFlags update_flags); /** * Constructor. Uses MappingQ1 * implicitly. */ - FEFaceValues (const FiniteElement& fe, - const Quadrature& quadrature, - const UpdateFlags update_flags); + FEFaceValues (const FiniteElement &fe, + const Quadrature &quadrature, + const UpdateFlags update_flags); /** * Reinitialize the gradients, Jacobi @@ -1158,9 +1718,44 @@ class FEFaceValues : public FEFaceValuesBase * number @p face_no of @p cell * and the given finite element. */ - void reinit (const typename DoFHandler::cell_iterator& cell, - const unsigned int face_no); - + void reinit (const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no); + + /** + * Reinitialize the gradients, + * Jacobi determinants, etc for + * the given face on given cell + * of type "iterator into a + * Triangulation object", and the + * given finite element. Since + * iterators into triangulation + * alone only convey information + * about the geometry of a cell, + * but not about degrees of + * freedom possibly associated + * with this cell, you will not + * be able to call some functions + * of this class if they need + * information about degrees of + * freedom. These functions are, + * above all, the + * get_function_value/grads/2nd_derivatives + * functions. If you want to call + * these functions, you have to + * call the @p reinit variants + * that take iterators into + * @p DoFHandler or other DoF + * handler type objects. + */ + void reinit (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no); + + /** + * Return a reference to this + * very object. + */ +//TODO: explain reason + const FEFaceValues & get_present_fe_values () const; private: /** @@ -1168,6 +1763,20 @@ class FEFaceValues : public FEFaceValuesBase * constructors. */ void initialize (const UpdateFlags update_flags); + + /** + * The reinit() functions do + * only that part of the work + * that requires knowledge of the + * type of iterator. After + * setting present_cell(), + * they pass on to this function, + * which does the real work, and + * which is independent of the + * actual type of the cell + * iterator. + */ + void do_reinit (const unsigned int face_no); }; @@ -1199,30 +1808,74 @@ class FESubfaceValues : public FEFaceValuesBase * matching the quadrature rule * and update flags. */ - FESubfaceValues (const Mapping& mapping, - const FiniteElement& fe, - const Quadrature& face_quadrature, - const UpdateFlags update_flags); + FESubfaceValues (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature &face_quadrature, + const UpdateFlags update_flags); /** * Constructor. Uses MappingQ1 * implicitly. */ - FESubfaceValues (const FiniteElement& fe, - const Quadrature& face_quadrature, - const UpdateFlags update_flags); - - /** - * Reinitialize the gradients, Jacobi - * determinants, etc for the face with - * number @p face_no of @p cell - * and the given finite element. - */ - void reinit (const typename DoFHandler::cell_iterator& cell, - const unsigned int face_no, - const unsigned int subface_no); + FESubfaceValues (const FiniteElement &fe, + const Quadrature &face_quadrature, + const UpdateFlags update_flags); + + /** + * Reinitialize the gradients, + * Jacobi determinants, etc for + * the given cell of type + * "iterator into a DoFHandler + * object", and the finite + * element associated with this + * object. It is assumed that the + * finite element used by the + * given cell is also the one + * used by this + * @p FESubfaceValues object. + */ + void reinit (const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no); + + /** + * Reinitialize the gradients, + * Jacobi determinants, etc for + * the given subface on given + * cell of type "iterator into a + * Triangulation object", and the + * given finite element. Since + * iterators into triangulation + * alone only convey information + * about the geometry of a cell, + * but not about degrees of + * freedom possibly associated + * with this cell, you will not + * be able to call some functions + * of this class if they need + * information about degrees of + * freedom. These functions are, + * above all, the + * get_function_value/grads/2nd_derivatives + * functions. If you want to call + * these functions, you have to + * call the @p reinit variants + * that take iterators into + * @p DoFHandler or other DoF + * handler type objects. + */ + void reinit (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no); + + /** + * Return a reference to this + * very object. + */ +//TODO: explain reason + const FESubfaceValues & get_present_fe_values () const; - /** + /** * Exception */ DeclException0 (ExcReinitCalledWithBoundaryFace); @@ -1239,6 +1892,21 @@ class FESubfaceValues : public FEFaceValuesBase * constructors. */ void initialize (const UpdateFlags update_flags); + + /** + * The reinit() functions do + * only that part of the work + * that requires knowledge of the + * type of iterator. After + * setting present_cell(), + * they pass on to this function, + * which does the real work, and + * which is independent of the + * actual type of the cell + * iterator. + */ + void do_reinit (const unsigned int face_no, + const unsigned int subface_no); }; /*@}*/ @@ -1551,10 +2219,10 @@ FEValuesBase::get_mapping () const template inline -typename Triangulation::cell_iterator +const typename Triangulation::cell_iterator FEValuesBase::get_cell () const { - return present_cell; + return *present_cell; } /*------------------------ Inline functions: FEValues ----------------------------*/ @@ -1569,6 +2237,16 @@ FEValues::get_quadrature () const } + +template +inline +const FEValues & +FEValues::get_present_fe_values () const +{ + return *this; +} + + /*------------------------ Inline functions: FEFaceValuesBase --------------------*/ @@ -1586,6 +2264,39 @@ FEFaceValuesBase::normal_vector (const unsigned int i) const } + +/*------------------------ Inline functions: FE*FaceValues --------------------*/ + +template +inline +const Quadrature & +FEFaceValuesBase::get_quadrature () const +{ + return quadrature; +} + + + +template +inline +const FEFaceValues & +FEFaceValues::get_present_fe_values () const +{ + return *this; +} + + + +template +inline +const FESubfaceValues & +FESubfaceValues::get_present_fe_values () const +{ + return *this; +} + + + template inline const Tensor<1,dim> & @@ -1601,14 +2312,5 @@ FEFaceValuesBase::boundary_form (const unsigned int i) const -template -inline -const Quadrature & -FEFaceValuesBase::get_quadrature () const -{ - return quadrature; -} - - #endif diff --git a/deal.II/deal.II/include/fe/hp_fe_values.h b/deal.II/deal.II/include/fe/hp_fe_values.h new file mode 100644 index 0000000000..09c7cb1bbc --- /dev/null +++ b/deal.II/deal.II/include/fe/hp_fe_values.h @@ -0,0 +1,455 @@ +//---------------------------- hp_fe_values.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003, 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- hp_fe_values.h --------------------------- +#ifndef __deal2__hp_fe_values_h +#define __deal2__hp_fe_values_h + +#include +#include +#include +#include + +#include +#include + +template class MappingQ1; +template class FiniteElement; + + +namespace internal +{ +/** + * 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. + */ + 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 = 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); + + 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 >, + 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 hpFE*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 hpFEValues it is equal to + * the first one, while for face integration it is one less. + * + * @author Wolfgang Bangerth, 2003 + */ + template + class hpFEValuesBase + { + public: + /** + * Constructor. Set the fields + * of this class to the values + * indicated by the parameters + * to the constructor. + */ + hpFEValuesBase (const Mapping &mapping, + const Quadrature &quadrature, + const UpdateFlags update_flags); + + /** + * Constructor. Set the fields + * of this class to the values + * indicated by the parameters + * to the constructor, and + * choose a @p MappingQ1 + * object for the mapping + * object. + */ + hpFEValuesBase (const Quadrature &quadrature, + const UpdateFlags update_flags); + + protected: + /** + * Pointer to a mapping object as + * provided to the constructor, + * or a default @p MappingQ1 + * object if none was given. + */ + const SmartPointer > mapping; + + /** + * Pointer to a quadrature object + * as provided to the + * constructor. + */ + const SmartPointer > quadrature; + + /** + * Values of the update flags as + * given to the constructor. + */ + const UpdateFlags update_flags; + + private: + /** + * Default mapping, in case + * none was provided through + * the call to the constructor. + */ + static const MappingQ1 default_mapping; + }; +} + + + +template +class hpFEValues : public internal::FEValuesMap >, + protected internal::hpFEValuesBase +{ + public: + /** + * Constructor. Initialize this + * object with the given + * parameters. + * + * The finite element + * collection parameter is + * actually ignored, but is in + * the signature of this + * function to make it + * compatible with the + * signature of the respective + * constructor of the usual + * FEValues object, with + * the respective parameter in + * that function also being the + * return value of the + * DoFHandler::get_fe() + * function. + */ + hpFEValues (const Mapping &mapping, + const FECollection &fe_collection, + const Quadrature &quadrature, + const UpdateFlags update_flags); + + /** + * Constructor. Initialize this + * object with the given + * parameters, and choose a + * @p MappingQ1 object for the + * mapping object. + * + * The finite element + * collection parameter is + * actually ignored, but is in + * the signature of this + * function to make it + * compatible with the + * signature of the respective + * constructor of the usual + * FEValues object, with + * the respective parameter in + * that function also being the + * return value of the + * DoFHandler::get_fe() + * function. + */ + hpFEValues (const FECollection &fe_collection, + const Quadrature &quadrature, + 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 hpDoFHandler::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; +}; + + + +template +class hpFEFaceValues : public internal::FEValuesMap >, + protected internal::hpFEValuesBase +{ + public: + /** + * Constructor. Initialize this + * object with the given + * parameters. + * + * The finite element + * collection parameter is + * actually ignored, but is in + * the signature of this + * function to make it + * compatible with the + * signature of the respective + * constructor of the usual + * FEValues object, with + * the respective parameter in + * that function also being the + * return value of the + * DoFHandler::get_fe() + * function. + */ + hpFEFaceValues (const Mapping &mapping, + const FECollection &fe_collection, + const Quadrature &quadrature, + const UpdateFlags update_flags); + + /** + * Constructor. Initialize this + * object with the given + * parameters, and choose a + * @p MappingQ1 object for the + * mapping object. + * + * The finite element + * collection parameter is + * actually ignored, but is in + * the signature of this + * function to make it + * compatible with the + * signature of the respective + * constructor of the usual + * FEValues object, with + * the respective parameter in + * that function also being the + * return value of the + * DoFHandler::get_fe() + * function. + */ + hpFEFaceValues (const FECollection &fe_collection, + const Quadrature &quadrature, + 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 hpDoFHandler::cell_iterator &cell, + const unsigned int face_no); + + protected: + /** + * Create an object of type + * @p FEValues for this + * particular finite element. + */ + virtual + FEFaceValues * + create_fe_values (const FiniteElement &fe) const; +}; + + + +template +class hpFESubfaceValues : public internal::FEValuesMap >, + protected internal::hpFEValuesBase +{ + public: + /** + * Constructor. Initialize this + * object with the given + * parameters. + * + * The finite element + * collection parameter is + * actually ignored, but is in + * the signature of this + * function to make it + * compatible with the + * signature of the respective + * constructor of the usual + * FEValues object, with + * the respective parameter in + * that function also being the + * return value of the + * DoFHandler::get_fe() + * function. + */ + hpFESubfaceValues (const Mapping &mapping, + const FECollection &fe_collection, + const Quadrature &quadrature, + const UpdateFlags update_flags); + + /** + * Constructor. Initialize this + * object with the given + * parameters, and choose a + * @p MappingQ1 object for the + * mapping object. + * + * The finite element + * collection parameter is + * actually ignored, but is in + * the signature of this + * function to make it + * compatible with the + * signature of the respective + * constructor of the usual + * FEValues object, with + * the respective parameter in + * that function also being the + * return value of the + * DoFHandler::get_fe() + * function. + */ + hpFESubfaceValues (const FECollection &fe_collection, + const Quadrature &quadrature, + 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 hpDoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no); + + protected: + /** + * Create an object of type + * @p FEValues for this + * particular finite element. + */ + virtual + FESubfaceValues * + create_fe_values (const FiniteElement &fe) const; +}; + + +// -------------- inline and template functions -------------- + +namespace internal +{ + template + inline + const FEValues & + FEValuesMap::get_present_fe_values () const + { + return *present_fe_values; + } +} + + + + +#endif diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index 5613229587..dcfbb476f1 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -15,13 +15,12 @@ #include -#include -#include -#include #include #include #include +#include + template class Quadrature; template class FEValuesData; template class FEValues; @@ -445,7 +444,7 @@ class Mapping : public Subscriptor * */ virtual void - fill_fe_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_values (const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, InternalDataBase &internal, std::vector > &quadrature_points, @@ -473,7 +472,7 @@ class Mapping : public Subscriptor * transformed quadrature weight. */ virtual void - fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_face_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &quadrature, InternalDataBase &internal, @@ -486,7 +485,7 @@ class Mapping : public Subscriptor * See above. */ virtual void - fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h index b174e0969f..0ad66dce35 100644 --- a/deal.II/deal.II/include/fe/mapping_cartesian.h +++ b/deal.II/deal.II/include/fe/mapping_cartesian.h @@ -54,14 +54,14 @@ class MappingCartesian : public Mapping const Quadrature& quadrature) const; virtual void - fill_fe_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_values (const typename Triangulation::cell_iterator &cell, const Quadrature& quadrature, typename Mapping::InternalDataBase &mapping_data, std::vector > &quadrature_points, std::vector &JxW_values) const ; virtual void - fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_face_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature& quadrature, typename Mapping::InternalDataBase &mapping_data, @@ -70,7 +70,7 @@ class MappingCartesian : public Mapping std::vector > &boundary_form, std::vector > &normal_vectors) const ; virtual void - fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature& quadrature, @@ -193,7 +193,7 @@ class MappingCartesian : public Mapping * Do the computation for the * fill_* functions. */ - void compute_fill (const typename DoFHandler::cell_iterator &cell, + void compute_fill (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, InternalData& data, @@ -216,7 +216,7 @@ class MappingCartesian : public Mapping /// @if NoDoc template <> void MappingCartesian<1>::fill_fe_face_values ( - const DoFHandler<1>::cell_iterator &, + const Triangulation<1>::cell_iterator &, const unsigned, const Quadrature<0>&, Mapping<1>::InternalDataBase&, @@ -226,7 +226,7 @@ template <> void MappingCartesian<1>::fill_fe_face_values ( std::vector >&) const; template <> void MappingCartesian<1>::fill_fe_subface_values ( - const DoFHandler<1>::cell_iterator &, + const Triangulation<1>::cell_iterator &, const unsigned, const unsigned, const Quadrature<0>&, diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index 100b4fc952..16b7d79d9f 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -184,7 +184,7 @@ class MappingQ : public MappingQ1 * Mapping. */ virtual void - fill_fe_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_values (const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, typename std::vector > &quadrature_points, @@ -195,7 +195,7 @@ class MappingQ : public MappingQ1 * Mapping. */ virtual void - fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_face_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature& quadrature, typename Mapping::InternalDataBase &mapping_data, @@ -209,7 +209,7 @@ class MappingQ : public MappingQ1 * Mapping. */ virtual void - fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature& quadrature, diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index c67cd38f12..3dfe781fbc 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -324,7 +324,7 @@ class MappingQ1 : public Mapping * @p mapping_support_points are * stored. */ - typename DoFHandler::cell_iterator cell_of_current_support_points; + typename Triangulation::cell_iterator cell_of_current_support_points; /** * Default value of this flag @@ -374,7 +374,7 @@ class MappingQ1 : public Mapping * Mapping. */ virtual void - fill_fe_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_values (const typename Triangulation::cell_iterator &cell, const Quadrature& quadrature, typename Mapping::InternalDataBase &mapping_data, typename std::vector > &quadrature_points, @@ -385,7 +385,7 @@ class MappingQ1 : public Mapping * Mapping. */ virtual void - fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_face_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature& quadrature, typename Mapping::InternalDataBase &mapping_data, @@ -399,7 +399,7 @@ class MappingQ1 : public Mapping * Mapping. */ virtual void - fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, + fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature& quadrature, @@ -458,7 +458,7 @@ class MappingQ1 : public Mapping * Do the computation for the * fill_* functions. */ - void compute_fill (const typename DoFHandler::cell_iterator &cell, + void compute_fill (const typename Triangulation::cell_iterator &cell, const unsigned int npts, const DataSetDescriptor data_set, InternalData &data, @@ -468,7 +468,7 @@ class MappingQ1 : public Mapping * Do the computation for the * fill_* functions. */ - void compute_fill_face (const typename DoFHandler::cell_iterator &cell, + void compute_fill_face (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const bool is_subface, const unsigned int npts, @@ -697,7 +697,7 @@ MappingQ1<1>::compute_face_data (const UpdateFlags, InternalData &) const; template <> void MappingQ1<1>::compute_fill_face ( - const DoFHandler<1>::cell_iterator &, + const Triangulation<1>::cell_iterator &, const unsigned int, const bool, const unsigned int, @@ -710,7 +710,7 @@ template <> void MappingQ1<1>::compute_fill_face ( std::vector > &) const; template <> void MappingQ1<1>::fill_fe_face_values ( - const DoFHandler<1>::cell_iterator &, + const Triangulation<1>::cell_iterator &, const unsigned, const Quadrature<0>&, Mapping<1>::InternalDataBase&, @@ -720,7 +720,7 @@ template <> void MappingQ1<1>::fill_fe_face_values ( std::vector >&) const; template <> void MappingQ1<1>::fill_fe_subface_values ( - const DoFHandler<1>::cell_iterator &, + const Triangulation<1>::cell_iterator &, const unsigned, const unsigned, const Quadrature<0>&, diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 93b87e4aab..0a83912056 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -12,9 +12,9 @@ //---------------------------- fe.cc --------------------------- -#include #include #include +#include #include #include #include @@ -462,7 +462,7 @@ template void FiniteElementBase:: compute_2nd (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int, typename Mapping::InternalDataBase &mapping_internal, InternalDataBase &fe_internal, diff --git a/deal.II/deal.II/source/fe/fe_data.cc b/deal.II/deal.II/source/fe/fe_data.cc index 5769609076..5d3f987a5b 100644 --- a/deal.II/deal.II/source/fe/fe_data.cc +++ b/deal.II/deal.II/source/fe/fe_data.cc @@ -12,9 +12,8 @@ //---------------------------------------------------------------- -#include #include - +#include template diff --git a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc index 32ba0c7bbb..e9111b188c 100644 --- a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc +++ b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc @@ -286,7 +286,7 @@ template void FE_DGPNonparametric::fill_fe_values ( const Mapping&, - const typename DoFHandler::cell_iterator&, + const typename Triangulation::cell_iterator&, const Quadrature&, typename Mapping::InternalDataBase&, typename Mapping::InternalDataBase& fedata, @@ -326,7 +326,7 @@ template void FE_DGPNonparametric::fill_fe_face_values ( const Mapping&, - const typename DoFHandler::cell_iterator&, + const typename Triangulation::cell_iterator&, const unsigned int, const Quadrature&, typename Mapping::InternalDataBase&, @@ -367,7 +367,7 @@ template void FE_DGPNonparametric::fill_fe_subface_values ( const Mapping&, - const typename DoFHandler::cell_iterator&, + const typename Triangulation::cell_iterator&, const unsigned int, const unsigned int, const Quadrature&, diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index c679cc49bf..bd7dc5fc4c 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -537,7 +537,6 @@ get_interpolation_matrix (const FiniteElementBase &x_source_fe, // Fill data of FEValues //---------------------------------------------------------------------- - template bool FE_DGQ::has_support_on_face (const unsigned int shape_index, diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index 4ae45426a6..92c7788f45 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -979,7 +979,7 @@ FE_Nedelec::get_data (const UpdateFlags update_flags, template void FE_Nedelec::fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fedata, @@ -1104,7 +1104,7 @@ FE_Nedelec::fill_fe_values (const Mapping &mapping, template void FE_Nedelec::fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, @@ -1249,7 +1249,7 @@ FE_Nedelec::fill_fe_face_values (const Mapping &mapp template void FE_Nedelec::fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face, const unsigned int subface, const Quadrature &quadrature, diff --git a/deal.II/deal.II/source/fe/fe_poly.cc b/deal.II/deal.II/source/fe/fe_poly.cc index 83e55ce5fa..38e8f75c1a 100644 --- a/deal.II/deal.II/source/fe/fe_poly.cc +++ b/deal.II/deal.II/source/fe/fe_poly.cc @@ -235,11 +235,11 @@ FE_Poly::get_data (const UpdateFlags update_flags, template void FE_Poly::fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename Mapping::InternalDataBase &fedata, - FEValuesData &data) const + const typename Triangulation::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { // convert data object to internal // data for this class. fails with @@ -277,7 +277,7 @@ FE_Poly::fill_fe_values (const Mapping &mapping template void FE_Poly::fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, @@ -327,7 +327,7 @@ FE_Poly::fill_fe_face_values (const Mapping &ma template void FE_Poly::fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face, const unsigned int subface, const Quadrature &quadrature, diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index d925e7cde3..8c820606b6 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -930,7 +930,6 @@ FE_Q<1>::initialize_constraints () #endif - #if deal_II_dimension == 2 template <> diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc index e045278fe1..1bc800305a 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -1493,7 +1493,7 @@ FE_RaviartThomas::get_data (const UpdateFlags update_flags, template void FE_RaviartThomas::fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fedata, @@ -1618,7 +1618,7 @@ FE_RaviartThomas::fill_fe_values (const Mapping &map template void FE_RaviartThomas::fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, @@ -1755,7 +1755,7 @@ FE_RaviartThomas::fill_fe_face_values (const Mapping template void FE_RaviartThomas::fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face, const unsigned int subface, const Quadrature &quadrature, diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index ad1acc86d2..4a21e06a25 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -670,7 +670,7 @@ template void FESystem:: fill_fe_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fe_data, @@ -686,7 +686,7 @@ template void FESystem:: fill_fe_face_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, @@ -704,7 +704,7 @@ template void FESystem:: fill_fe_subface_values (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, @@ -723,7 +723,7 @@ template void FESystem:: compute_fill (const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, + const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &quadrature, diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 28eeb32df8..c6117959e6 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -34,6 +34,8 @@ #include #include +#include + #ifdef HAVE_STD_STRINGSTREAM # include #else @@ -44,6 +46,30 @@ namespace { + template + unsigned int + max_dofs_per_cell (const DoFHandler &dh) + { + return dh.get_fe().dofs_per_cell; + } + + + template + unsigned int + max_dofs_per_face (const DoFHandler &dh) + { + return dh.get_fe().dofs_per_face; + } + + + template + unsigned int + n_components (const DoFHandler &dh) + { + return dh.get_fe().n_components(); + } + + // forwarder function for // FE::get_interpolation_matrix. we // will want to call that function @@ -76,8 +102,6 @@ namespace interpolation_matrix = tmp; } - - // return true if the given pattern // string matches the given name at // the first position of the string @@ -352,6 +376,7 @@ void FETools::get_interpolation_difference_matrix (const FiniteElement &fe1 } + template void FETools::get_projection_matrix (const FiniteElement &fe1, const FiniteElement &fe2, @@ -444,11 +469,16 @@ void FETools::get_projection_matrix (const FiniteElement &fe1, } -template -void FETools::interpolate(const DoFHandler &dof1, - const InVector &u1, - const DoFHandler &dof2, - OutVector &u2) + +template class DH1, + template class DH2, + class InVector, class OutVector> +void +FETools::interpolate(const DH1 &dof1, + const InVector &u1, + const DH2 &dof2, + OutVector &u2) { ConstraintMatrix dummy; dummy.close(); @@ -457,17 +487,21 @@ void FETools::interpolate(const DoFHandler &dof1, template -void FETools::interpolate(const DoFHandler &dof1, - const InVector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints, - OutVector &u2) +void +FETools::interpolate(const DoFHandler &dof1, + const InVector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints, + OutVector &u2) { - Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch()); + Assert(&dof1.get_tria() == &dof2.get_tria(), ExcTriangulationMismatch()); Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); - Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs())); + ExcDimensionMismatch(dof1.get_fe().n_components(), + dof2.get_fe().n_components())); + Assert(u1.size()==dof1.n_dofs(), + ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u2.size()==dof2.n_dofs(), + ExcDimensionMismatch(u2.size(), dof2.n_dofs())); // for continuous elements on grids // with hanging nodes we need @@ -475,9 +509,9 @@ void FETools::interpolate(const DoFHandler &dof1, // constraints. Consequentely, if // there are no constraints then // hanging nodes are not allowed. - const bool hanging_nodes_not_allowed= - (dof2.get_fe().dofs_per_vertex != 0) - && (constraints.n_constraints()==0); + const bool hanging_nodes_not_allowed + = ((dof2.get_fe().dofs_per_vertex != 0) && + (constraints.n_constraints() == 0)); const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell; const unsigned int dofs_per_cell2=dof2.get_fe().dofs_per_cell; @@ -495,11 +529,11 @@ void FETools::interpolate(const DoFHandler &dof1, cell2 = dof2.begin_active(), endc2 = dof2.end(); - std::vector index_multiplicity(dof2.n_dofs(),0); + std::vector touch_count(dof2.n_dofs(),0); std::vector dofs (dofs_per_cell2); u2 = 0; - for (; cell1!=endc1, cell2!=endc2; ++cell1, ++cell2) + for (; cell1!=endc1; ++cell1, ++cell2) { if (hanging_nodes_not_allowed) for (unsigned int face=0; face::faces_per_cell; ++face) @@ -513,17 +547,142 @@ void FETools::interpolate(const DoFHandler &dof1, for (unsigned int i=0; i class DH1, + template class DH2, + class InVector, class OutVector> +void +FETools::interpolate (const DH1 &dof1, + const InVector &u1, + const DH2 &dof2, + const ConstraintMatrix &constraints, + OutVector &u2) +{ + Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch()); + + Assert(u1.size()==dof1.n_dofs(), + ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u2.size()==dof2.n_dofs(), + ExcDimensionMismatch(u2.size(), dof2.n_dofs())); + + // allocate vectors at maximal + // size. will be reinited in inner + // cell, but Vector makes sure that + // this does not lead to + // reallocation of memory + Vector u1_local(max_dofs_per_cell(dof1)); + Vector u2_local(max_dofs_per_cell(dof2)); + + // have a map for interpolation + // matrices. shared_ptr make sure + // that memory is released again + std::map *, + std::map *, + boost::shared_ptr > > > + interpolation_matrices; + + typename DH1::active_cell_iterator cell1 = dof1.begin_active(), + endc1 = dof1.end(); + typename DH2::active_cell_iterator cell2 = dof2.begin_active(), + endc2 = dof2.end(); + + std::vector touch_count(dof2.n_dofs(),0); + std::vector dofs; + dofs.reserve (max_dofs_per_cell (dof2)); + u2.clear (); + + for (; cell1!=endc1; ++cell1, ++cell2) + { + Assert(cell1->get_fe().n_components() == cell2->get_fe().n_components(), + ExcDimensionMismatch (cell1->get_fe().n_components(), + cell2->get_fe().n_components())); + + // for continuous elements on + // grids with hanging nodes we + // need hanging node + // constraints. Consequentely, + // if there are no constraints + // then hanging nodes are not + // allowed. + const bool hanging_nodes_not_allowed + = ((cell2->get_fe().dofs_per_vertex != 0) && + (constraints.n_constraints() == 0)); + + if (hanging_nodes_not_allowed) + for (unsigned int face=0; face::faces_per_cell; ++face) + Assert (cell1->at_boundary(face) || + cell1->neighbor(face)->level() == cell1->level(), + ExcHangingNodesNotAllowed(0)); + + + const unsigned int dofs_per_cell1 = cell1->get_fe().dofs_per_cell; + const unsigned int dofs_per_cell2 = cell2->get_fe().dofs_per_cell; + u1_local.reinit (dofs_per_cell1); + u2_local.reinit (dofs_per_cell2); + + // check if interpolation + // matrix for this particular + // pair of elements is already + // there + if (interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()] == 0) + { + boost::shared_ptr > + interpolation_matrix (new FullMatrix (dofs_per_cell2, + dofs_per_cell1)); + interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()] + = interpolation_matrix; + + FETools::get_interpolation_matrix(cell1->get_fe(), + cell2->get_fe(), + *interpolation_matrix); + } + + cell1->get_dof_values(u1, u1_local); + interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()] + ->vmult(u2_local, u1_local); + + dofs.resize (dofs_per_cell2); + cell2->get_dof_indices(dofs); + + for (unsigned int i=0; i &dof1, } + template -void FETools::back_interpolate(const DoFHandler &dof1, - const InVector &u1, - const FiniteElement &fe2, - OutVector &u1_interpolated) +void +FETools::back_interpolate(const DoFHandler &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_interpolated) { Assert(dof1.get_fe().n_components() == fe2.n_components(), ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components())); @@ -545,7 +706,7 @@ void FETools::back_interpolate(const DoFHandler &dof1, // For continuous elements on grids // with hanging nodes we need - // hnaging node + // hanging node // constraints. Consequently, when // the elements are continuous no // hanging node constraints are @@ -580,6 +741,78 @@ void FETools::back_interpolate(const DoFHandler &dof1, +template class DH, + class InVector, class OutVector> +void +FETools::back_interpolate(const DH &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_interpolated) +{ + Assert(u1.size() == dof1.n_dofs(), + ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u1_interpolated.size() == dof1.n_dofs(), + ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); + + Vector u1_local(max_dofs_per_cell(dof1)); + Vector u1_int_local(max_dofs_per_cell(dof1)); + + typename DH::active_cell_iterator cell = dof1.begin_active(), + endc = dof1.end(); + + // map from possible fe objects in + // dof1 to the back_interpolation + // matrices + std::map *, + boost::shared_ptr > > interpolation_matrices; + + for (; cell!=endc; ++cell) + { + Assert(cell->get_fe().n_components() == fe2.n_components(), + ExcDimensionMismatch(cell->get_fe().n_components(), + fe2.n_components())); + + // For continuous elements on + // grids with hanging nodes we + // need hanging node + // constraints. Consequently, + // when the elements are + // continuous no hanging node + // constraints are allowed. + const bool hanging_nodes_not_allowed= + (cell->get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); + + if (hanging_nodes_not_allowed) + for (unsigned int face=0; face::faces_per_cell; ++face) + Assert (cell->at_boundary(face) || + cell->neighbor(face)->level() == cell->level(), + ExcHangingNodesNotAllowed(0)); + + const unsigned int dofs_per_cell1 = cell->get_fe().dofs_per_cell; + + // make sure back_interpolation + // matrix is available + if (interpolation_matrices[&cell->get_fe()] != 0) + { + interpolation_matrices[&cell->get_fe()] = + boost::shared_ptr > + (new FullMatrix(dofs_per_cell1, dofs_per_cell1)); + get_back_interpolation_matrix(dof1.get_fe(), fe2, + *interpolation_matrices[&cell->get_fe()]); + } + + u1_local.reinit (dofs_per_cell1); + u1_int_local.reinit (dofs_per_cell1); + + cell->get_dof_values(u1, u1_local); + interpolation_matrices[&cell->get_fe()]->vmult(u1_int_local, u1_local); + cell->set_dof_values(u1_int_local, u1_interpolated); + } +} + + + template void FETools::back_interpolate(const DoFHandler &dof1, const ConstraintMatrix &constraints1, diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index bc65ab2486..6ae075c98e 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -12,9 +12,6 @@ //---------------------------- fe_values.cc --------------------------- -#include -#include -#include #include #include #include @@ -25,10 +22,241 @@ #include #include #include +#include +#include +#include #include +/* ---------------- FEValuesBase::CellIteratorBase --------- */ + +template +inline +FEValuesBase::CellIteratorBase::~CellIteratorBase () +{} + + + +/* ---------------- FEValuesBase::CellIterator --------- */ + + +template +template +inline +FEValuesBase::CellIterator::CellIterator (const CI &cell) + : + cell(cell) +{} + + + +template +template +FEValuesBase::CellIterator:: +operator const typename Triangulation::cell_iterator () const +{ + return cell; +} + + + +template +template +unsigned int +FEValuesBase::CellIterator::n_dofs_for_dof_handler () const +{ + return cell->get_dof_handler().n_dofs(); +} + + + +template +template +void +FEValuesBase::CellIterator:: +get_interpolated_dof_values (const Vector &in, + Vector &out) const +{ + cell->get_interpolated_dof_values (in, out); +} + + +template +template +void +FEValuesBase::CellIterator:: +get_interpolated_dof_values (const Vector &in, + Vector &out) const +{ + cell->get_interpolated_dof_values (in, out); +} + + +template +template +void +FEValuesBase::CellIterator:: +get_interpolated_dof_values (const BlockVector &in, + Vector &out) const +{ + cell->get_interpolated_dof_values (in, out); +} + + +template +template +void +FEValuesBase::CellIterator:: +get_interpolated_dof_values (const BlockVector &in, + Vector &out) const +{ + cell->get_interpolated_dof_values (in, out); +} + + + +#ifdef DEAL_II_USE_PETSC + +template +template +void +FEValuesBase::CellIterator:: +get_interpolated_dof_values (const PETScWrappers::Vector &in, + Vector &out) const +{ + cell->get_interpolated_dof_values (in, out); +} + + + +template +template +void +FEValuesBase::CellIterator:: +get_interpolated_dof_values (const PETScWrappers::BlockVector &in, + Vector &out) const +{ + cell->get_interpolated_dof_values (in, out); +} + +#endif + + +/* ---------------- FEValuesBase::TriaCellIterator --------- */ + +template +const char * const +FEValuesBase::TriaCellIterator::message_string += ("You have previously called the FEValues::reinit function with a\n" + "cell iterator of type Triangulation::cell_iterator. However,\n" + "when you do this, you cannot call some functions in the FEValues\n" + "class, such as the get_function_values/grads/2nd_derivatives\n" + "functions. If you need these functions, then you need to call\n" + "FEValues::reinit with an iterator type that allows to extract\n" + "degrees of freedom, such as DoFHandler::cell_iterator."); + + +template +inline +FEValuesBase::TriaCellIterator:: +TriaCellIterator (const typename Triangulation::cell_iterator &cell) + : + cell(cell) +{} + + + +template +FEValuesBase::TriaCellIterator:: +operator const typename Triangulation::cell_iterator () const +{ + return cell; +} + + + +template +unsigned int +FEValuesBase::TriaCellIterator::n_dofs_for_dof_handler () const +{ + Assert (false, ExcMessage (message_string)); + return 0; +} + + +template +void +FEValuesBase::TriaCellIterator:: +get_interpolated_dof_values (const Vector &, + Vector &) const +{ + Assert (false, ExcMessage (message_string)); +} + + +template +void +FEValuesBase::TriaCellIterator:: +get_interpolated_dof_values (const Vector &, + Vector &) const +{ + Assert (false, ExcMessage (message_string)); +} + + +template +void +FEValuesBase::TriaCellIterator:: +get_interpolated_dof_values (const BlockVector &, + Vector &) const +{ + Assert (false, ExcMessage (message_string)); +} + + + +template +void +FEValuesBase::TriaCellIterator:: +get_interpolated_dof_values (const BlockVector &, + Vector &) const +{ + Assert (false, ExcMessage (message_string)); +} + + + +#ifdef DEAL_II_USE_PETSC + +template +void +FEValuesBase::TriaCellIterator:: +get_interpolated_dof_values (const PETScWrappers::Vector &, + Vector &) const +{ + Assert (false, ExcMessage (message_string)); +} + + + +template +void +FEValuesBase::TriaCellIterator:: +get_interpolated_dof_values (const PETScWrappers::BlockVector &, + Vector &) const +{ + Assert (false, ExcMessage (message_string)); +} + +#endif + + + + +/* --------------------- FEValuesData ----------------- */ + + template void FEValuesData::initialize (const unsigned int n_quadrature_points, @@ -144,16 +372,13 @@ void FEValuesBase::get_function_values (const InputVector &fe_fu ExcWrongNoOfComponents()); Assert (values.size() == n_quadrature_points, ExcWrongVectorSize(values.size(), n_quadrature_points)); - Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), - ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); + Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(), + ExcWrongVectorSize(fe_function.size(), present_cell->n_dofs_for_dof_handler())); // get function values of dofs // on this cell Vector dof_values (dofs_per_cell); - if (present_cell->active()) - present_cell->get_dof_values (fe_function, dof_values); - else - present_cell->get_interpolated_dof_values(fe_function, dof_values); + present_cell->get_interpolated_dof_values(fe_function, dof_values); // initialize with zero std::fill_n (values.begin(), n_quadrature_points, 0); @@ -185,16 +410,13 @@ void FEValuesBase::get_function_values (const InputVector Assert (values[i].size() == n_components, ExcWrongNoOfComponents()); Assert (this->update_flags & update_values, ExcAccessToUninitializedField()); - Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), - ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); + Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(), + ExcWrongVectorSize(fe_function.size(), present_cell->n_dofs_for_dof_handler())); // get function values of dofs // on this cell Vector dof_values (dofs_per_cell); - if (present_cell->active()) - present_cell->get_dof_values (fe_function, dof_values); - else - present_cell->get_interpolated_dof_values(fe_function, dof_values); + present_cell->get_interpolated_dof_values(fe_function, dof_values); // initialize with zero for (unsigned i=0;iget_dof_handler().n_dofs(), - ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); + Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(), + ExcWrongVectorSize(fe_function.size(), present_cell->n_dofs_for_dof_handler())); // get function values of dofs // on this cell Vector dof_values (dofs_per_cell); - if (present_cell->active()) - present_cell->get_dof_values (fe_function, dof_values); - else - present_cell->get_interpolated_dof_values(fe_function, dof_values); + present_cell->get_interpolated_dof_values(fe_function, dof_values); // initialize with zero std::fill_n (gradients.begin(), n_quadrature_points, Tensor<1,dim>()); @@ -298,16 +517,13 @@ get_function_grads (const InputVector &fe_function, Assert (gradients[i].size() == n_components, ExcWrongNoOfComponents()); Assert (this->update_flags & update_gradients, ExcAccessToUninitializedField()); - Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), - ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); + Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(), + ExcWrongVectorSize(fe_function.size(), present_cell->n_dofs_for_dof_handler())); // get function values of dofs // on this cell Vector dof_values (dofs_per_cell); - if (present_cell->active()) - present_cell->get_dof_values (fe_function, dof_values); - else - present_cell->get_interpolated_dof_values(fe_function, dof_values); + present_cell->get_interpolated_dof_values(fe_function, dof_values); // initialize with zero for (unsigned i=0;iupdate_flags & update_second_derivatives, ExcAccessToUninitializedField()); - Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), - ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); + Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(), + ExcWrongVectorSize(fe_function.size(), present_cell->n_dofs_for_dof_handler())); // get function values of dofs // on this cell Vector dof_values (dofs_per_cell); - if (present_cell->active()) - present_cell->get_dof_values (fe_function, dof_values); - else - present_cell->get_interpolated_dof_values(fe_function, dof_values); + present_cell->get_interpolated_dof_values(fe_function, dof_values); // initialize with zero std::fill_n (second_derivatives.begin(), n_quadrature_points, Tensor<2,dim>()); @@ -397,16 +610,13 @@ get_function_2nd_derivatives (const InputVector &fe_func Assert (second_derivs[i].size() == n_components, ExcWrongNoOfComponents()); Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField()); - Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), - ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); + Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(), + ExcWrongVectorSize(fe_function.size(), present_cell->n_dofs_for_dof_handler())); // get function values of dofs // on this cell Vector dof_values (dofs_per_cell); - if (present_cell->active()) - present_cell->get_dof_values (fe_function, dof_values); - else - present_cell->get_interpolated_dof_values(fe_function, dof_values); + present_cell->get_interpolated_dof_values(fe_function, dof_values); // initialize with zero for (unsigned i=0;i::memory_consumption () const MemoryConsumption::memory_consumption (this->normal_vectors) + MemoryConsumption::memory_consumption (this->boundary_forms) + sizeof(this->update_flags) + - MemoryConsumption::memory_consumption (present_cell) + MemoryConsumption::memory_consumption (n_quadrature_points) + MemoryConsumption::memory_consumption (dofs_per_cell) + MemoryConsumption::memory_consumption (mapping) + @@ -580,26 +789,101 @@ FEValues::initialize (const UpdateFlags update_flags) template -void FEValues::reinit (const typename DoFHandler::cell_iterator &cell) +void +FEValues::reinit (const typename DoFHandler::cell_iterator &cell) { // assert that the finite elements // passed to the constructor and // used by the DoFHandler used by // this cell, are the same Assert (static_cast&>(*this->fe) == - static_cast&>(cell->get_dof_handler().get_fe()), + static_cast&>(cell->get_fe()), + typename FEValuesBase::ExcFEDontMatch()); + + // set new cell. auto_ptr will take + // care that old object gets + // destroyed and also that this + // object gets destroyed in the + // destruction of this class + this->present_cell.reset + (new typename FEValuesBase::template + CellIterator::cell_iterator> (cell)); + + // this was the part of the work + // that is dependent on the actual + // data type of the iterator. now + // pass on to the function doing + // the real work. + do_reinit (); +} + + + +template +void +FEValues::reinit (const typename MGDoFHandler::cell_iterator &cell) +{ + // assert that the finite elements + // passed to the constructor and + // used by the DoFHandler used by + // this cell, are the same + Assert (static_cast&>(*this->fe) == + static_cast&>(cell->get_fe()), typename FEValuesBase::ExcFEDontMatch()); - this->present_cell = cell; + // set new cell. auto_ptr will take + // care that old object gets + // destroyed and also that this + // object gets destroyed in the + // destruction of this class + this->present_cell.reset + (new typename FEValuesBase::template + CellIterator::cell_iterator> (cell)); + + // this was the part of the work + // that is dependent on the actual + // data type of the iterator. now + // pass on to the function doing + // the real work. + do_reinit (); +} + + + +template +void FEValues::reinit (const typename Triangulation::cell_iterator &cell) +{ + // no FE in this cell, so no check + // necessary here + + // set new cell. auto_ptr will take + // care that old object gets + // destroyed and also that this + // object gets destroyed in the + // destruction of this class + this->present_cell.reset + (new typename FEValuesBase::TriaCellIterator (cell)); + // this was the part of the work + // that is dependent on the actual + // data type of the iterator. now + // pass on to the function doing + // the real work. + do_reinit (); +} + + - this->get_mapping().fill_fe_values(cell, +template +void FEValues::do_reinit () +{ + this->get_mapping().fill_fe_values(*this->present_cell, quadrature, *this->mapping_data, this->quadrature_points, this->JxW_values); this->get_fe().fill_fe_values(this->get_mapping(), - cell, + *this->present_cell, quadrature, *this->mapping_data, *this->fe_data, @@ -740,9 +1024,57 @@ void FEFaceValues::reinit (const typename DoFHandler::cell_iterator &c static_cast&>(cell->get_dof_handler().get_fe()), typename FEValuesBase::ExcFEDontMatch()); - this->present_cell = cell; + Assert (face_no < GeometryInfo::faces_per_cell, + ExcIndexRange (face_no, 0, GeometryInfo::faces_per_cell)); - this->get_mapping().fill_fe_face_values(cell, face_no, + // set new cell. auto_ptr will take + // care that old object gets + // destroyed and also that this + // object gets destroyed in the + // destruction of this class + this->present_cell.reset + (new typename FEValuesBase::template + CellIterator::cell_iterator> (cell)); + + // this was the part of the work + // that is dependent on the actual + // data type of the iterator. now + // pass on to the function doing + // the real work. + do_reinit (face_no); +} + + + +template +void FEFaceValues::reinit (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no) +{ + Assert (face_no < GeometryInfo::faces_per_cell, + ExcIndexRange (face_no, 0, GeometryInfo::faces_per_cell)); + + // set new cell. auto_ptr will take + // care that old object gets + // destroyed and also that this + // object gets destroyed in the + // destruction of this class + this->present_cell.reset + (new typename FEValuesBase::TriaCellIterator (cell)); + + // this was the part of the work + // that is dependent on the actual + // data type of the iterator. now + // pass on to the function doing + // the real work. + do_reinit (face_no); +} + + + +template +void FEFaceValues::do_reinit (const unsigned int face_no) +{ + this->get_mapping().fill_fe_face_values(*this->present_cell, face_no, this->quadrature, *this->mapping_data, this->quadrature_points, @@ -751,7 +1083,7 @@ void FEFaceValues::reinit (const typename DoFHandler::cell_iterator &c this->normal_vectors); this->get_fe().fill_fe_face_values(this->get_mapping(), - cell, face_no, + *this->present_cell, face_no, this->quadrature, *this->mapping_data, *this->fe_data, @@ -841,9 +1173,59 @@ void FESubfaceValues::reinit (const typename DoFHandler::cell_iterator "already refined. Iterate over their children " "instead in these cases.")); - this->present_cell = cell; + // set new cell. auto_ptr will take + // care that old object gets + // destroyed and also that this + // object gets destroyed in the + // destruction of this class + this->present_cell.reset + (new typename FEValuesBase::template + CellIterator::cell_iterator> (cell)); + + // this was the part of the work + // that is dependent on the actual + // data type of the iterator. now + // pass on to the function doing + // the real work. + do_reinit (face_no, subface_no); +} + + - this->get_mapping().fill_fe_subface_values(cell, face_no, subface_no, +template +void FESubfaceValues::reinit (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no) +{ + Assert (face_no < GeometryInfo::faces_per_cell, + ExcIndexRange (face_no, 0, GeometryInfo::faces_per_cell)); + Assert (subface_no < GeometryInfo::subfaces_per_face, + ExcIndexRange (subface_no, 0, GeometryInfo::subfaces_per_face)); + + // set new cell. auto_ptr will take + // care that old object gets + // destroyed and also that this + // object gets destroyed in the + // destruction of this class + this->present_cell.reset + (new typename FEValuesBase::TriaCellIterator (cell)); + + // this was the part of the work + // that is dependent on the actual + // data type of the iterator. now + // pass on to the function doing + // the real work. + do_reinit (face_no, subface_no); +} + + + +template +void FESubfaceValues::do_reinit (const unsigned int face_no, + const unsigned int subface_no) +{ + this->get_mapping().fill_fe_subface_values(*this->present_cell, + face_no, subface_no, this->quadrature, *this->mapping_data, this->quadrature_points, @@ -852,7 +1234,8 @@ void FESubfaceValues::reinit (const typename DoFHandler::cell_iterator this->normal_vectors); this->get_fe().fill_fe_subface_values(this->get_mapping(), - cell, face_no, subface_no, + *this->present_cell, + face_no, subface_no, this->quadrature, *this->mapping_data, *this->fe_data, diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index 587cfa98ad..3c86eba344 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -127,7 +127,7 @@ MappingCartesian::get_subface_data (const UpdateFlags update_flags, template void -MappingCartesian::compute_fill (const typename DoFHandler::cell_iterator &cell, +MappingCartesian::compute_fill (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, InternalData &data, @@ -276,7 +276,7 @@ MappingCartesian::compute_fill (const typename DoFHandler::cell_iterat template void MappingCartesian:: -fill_fe_values (const typename DoFHandler::cell_iterator& cell, +fill_fe_values (const typename Triangulation::cell_iterator& cell, const Quadrature& q, typename Mapping::InternalDataBase& mapping_data, std::vector >& quadrature_points, @@ -314,7 +314,7 @@ fill_fe_values (const typename DoFHandler::cell_iterator& cell, template void -MappingCartesian::fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, +MappingCartesian::fill_fe_face_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &q, typename Mapping::InternalDataBase &mapping_data, @@ -356,7 +356,7 @@ MappingCartesian::fill_fe_face_values (const typename DoFHandler::cell template void -MappingCartesian::fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, +MappingCartesian::fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &q, @@ -400,7 +400,7 @@ MappingCartesian::fill_fe_subface_values (const typename DoFHandler::c template <> void -MappingCartesian<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &, +MappingCartesian<1>::fill_fe_face_values (const Triangulation<1>::cell_iterator &, const unsigned, const Quadrature<0>&, Mapping<1>::InternalDataBase&, @@ -415,7 +415,7 @@ MappingCartesian<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &, template <> void -MappingCartesian<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterator &, +MappingCartesian<1>::fill_fe_subface_values (const Triangulation<1>::cell_iterator &, const unsigned, const unsigned, const Quadrature<0>&, diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index 43456f694b..112196bd3c 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -258,7 +258,7 @@ MappingQ::get_subface_data (const UpdateFlags update_flags, template void -MappingQ::fill_fe_values (const typename DoFHandler::cell_iterator &cell, +MappingQ::fill_fe_values (const typename Triangulation::cell_iterator &cell, const Quadrature &q, typename Mapping::InternalDataBase &mapping_data, std::vector > &quadrature_points, @@ -295,7 +295,7 @@ MappingQ::fill_fe_values (const typename DoFHandler::cell_iterator &ce template void -MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, +MappingQ::fill_fe_face_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &q, typename Mapping::InternalDataBase &mapping_data, @@ -351,7 +351,7 @@ MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterato template void -MappingQ::fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, +MappingQ::fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &q, diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index 22abf10b9f..f9ceb4bb6d 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -439,7 +439,7 @@ MappingQ1::get_subface_data (const UpdateFlags update_flags, template void -MappingQ1::compute_fill (const typename DoFHandler::cell_iterator &cell, +MappingQ1::compute_fill (const typename Triangulation::cell_iterator &cell, const unsigned int n_q_points, const DataSetDescriptor data_set, InternalData &data, @@ -542,7 +542,7 @@ MappingQ1::compute_mapping_support_points( template void -MappingQ1::fill_fe_values (const typename DoFHandler::cell_iterator &cell, +MappingQ1::fill_fe_values (const typename Triangulation::cell_iterator &cell, const Quadrature &q, typename Mapping::InternalDataBase &mapping_data, std::vector > &quadrature_points, @@ -577,7 +577,7 @@ MappingQ1::fill_fe_values (const typename DoFHandler::cell_iterator &c template void -MappingQ1::compute_fill_face (const typename DoFHandler::cell_iterator &cell, +MappingQ1::compute_fill_face (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const bool is_subface, const unsigned int n_q_points, @@ -668,7 +668,7 @@ MappingQ1::compute_fill_face (const typename DoFHandler::cell_iterator template void -MappingQ1::fill_fe_face_values (const typename DoFHandler::cell_iterator &cell, +MappingQ1::fill_fe_face_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &q, typename Mapping::InternalDataBase &mapping_data, @@ -699,7 +699,7 @@ MappingQ1::fill_fe_face_values (const typename DoFHandler::cell_iterat template void -MappingQ1::fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, +MappingQ1::fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature &q, @@ -733,7 +733,7 @@ MappingQ1::fill_fe_subface_values (const typename DoFHandler::cell_ite template <> void -MappingQ1<1>::compute_fill_face (const DoFHandler<1>::cell_iterator &, +MappingQ1<1>::compute_fill_face (const Triangulation<1>::cell_iterator &, const unsigned int, const bool, const unsigned int, @@ -751,7 +751,7 @@ MappingQ1<1>::compute_fill_face (const DoFHandler<1>::cell_iterator &, template <> void -MappingQ1<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &, +MappingQ1<1>::fill_fe_face_values (const Triangulation<1>::cell_iterator &, const unsigned, const Quadrature<0>&, Mapping<1>::InternalDataBase&, @@ -766,7 +766,7 @@ MappingQ1<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &, template <> void -MappingQ1<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterator &, +MappingQ1<1>::fill_fe_subface_values (const Triangulation<1>::cell_iterator &, const unsigned, const unsigned, const Quadrature<0>&, diff --git a/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc b/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc index 40c83268a6..67c5a66373 100644 --- a/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc +++ b/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc @@ -23,7 +23,7 @@ template MappingQ1Eulerian::MappingQ1Eulerian (const Vector &euler_transform_vectors, - const DoFHandler &shiftmap_dof_handler) + const DoFHandler &shiftmap_dof_handler) : euler_transform_vectors(euler_transform_vectors), shiftmap_dof_handler(&shiftmap_dof_handler)