From: bangerth Date: Fri, 5 Dec 2008 22:27:44 +0000 (+0000) Subject: Merge branch_codimension_one. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fb428c78dfece7af8866fcfe651c0594bc640919;p=dealii-svn.git Merge branch_codimension_one. git-svn-id: https://svn.dealii.org/trunk@17866 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/geometry_info.h b/deal.II/base/include/base/geometry_info.h index 2296ee3df8..5bb4175122 100644 --- a/deal.II/base/include/base/geometry_info.h +++ b/deal.II/base/include/base/geometry_info.h @@ -810,6 +810,16 @@ struct GeometryInfo<0> */ static const unsigned int max_children_per_face = 0; + /** + * Return the number of children + * of a cell (or face) refined + * with ref_case. Since + * we are concerned here with + * points, the number of children + * is equal to one. + */ + static unsigned int n_children(const RefinementCase<0> &refinement_case); + /** * Number of vertices a cell has. */ @@ -1684,7 +1694,7 @@ struct GeometryInfo * entry. * * In the case of anisotropically refined - * cells and faces, the @p RefineCase of + * cells and faces, the @p RefinementCase of * the face, face_ref_case, * might have an influence on * which child is behind which given diff --git a/deal.II/base/include/base/quadrature_lib.h b/deal.II/base/include/base/quadrature_lib.h index 32f143b74e..0606d258ea 100644 --- a/deal.II/base/include/base/quadrature_lib.h +++ b/deal.II/base/include/base/quadrature_lib.h @@ -306,6 +306,72 @@ class QWeddle : public Quadrature QWeddle (); }; + + +/** + * Gauss Quadrature Formula to integrate ln[(x-a)/(b-a)]*f(x) on the + * interval [a,b], where f is a smooth function without + * singularities. See Numerical Recipes for more information. The + * quadrature formula weights are chosen in order to integrate + * f(x)*ln[(x-a)/(b-a)] (that is, the argument of the logarithm goes + * linearly from 0 to 1 on the interval of integration). So, the only + * thing one can specify is the function f(x). Using logarithm + * properties, one can add the integral of f(x)*ln(b-a) (calculated + * with QGauss class formulas) to the integral of f(x)*ln[(x-a)/(b-a)] + * to finally obtain the integral of f(x)*ln(x-a). It is also + * possible to integrate with weighting function ln[(b-x)/(b-a)] if + * revert is set to true. + * + * It works up to n=12. + * +//TODO: implement the calculation of points and weights as it is +//described in Numerical Recipes, instead of passing them to set_* +//functions, so that n can also be grater than 12. + * + */ +template +class QGaussLog : public Quadrature +{ + public: + /** + * Generate a formula with + * n quadrature points + */ + QGaussLog(const unsigned int n, + const bool revert=false); + + protected: + /** + * Sets the points of the + * quadrature formula. + */ + std::vector + set_quadrature_points(const unsigned int n) const; + + /** + * Sets the weights of the + * quadrature formula. + */ + std::vector + set_quadrature_weights(const unsigned int n) const; + + /** + * If the singularity is at the rightmost point + * of the interval just revert the vectors + * of points and weights and everything + * works. + */ + void + revert_points_and_weights(std::vector &p, + std::vector &w) const; + +}; + + + + + + /*@}*/ /* -------------- declaration of explicit specializations ------------- */ @@ -324,6 +390,11 @@ JacobiP(const long double, const int, const int, const unsigned int) const; template <> unsigned int QGaussLobatto<1>::gamma(const unsigned int n) const; + +template <> void QGaussLog<1>::revert_points_and_weights(std::vector &, std::vector &) const; +template <> std::vector QGaussLog<1>::set_quadrature_points(const unsigned int) const; +template <> std::vector QGaussLog<1>::set_quadrature_weights(const unsigned int) const; + template <> QGauss2<1>::QGauss2 (); template <> QGauss3<1>::QGauss3 (); template <> QGauss4<1>::QGauss4 (); @@ -335,6 +406,10 @@ template <> QTrapez<1>::QTrapez (); template <> QSimpson<1>::QSimpson (); template <> QMilne<1>::QMilne (); template <> QWeddle<1>::QWeddle (); +template <> QGaussLog<1>::QGaussLog (const unsigned int n, const bool revert); + + + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/base/include/base/tensor.h b/deal.II/base/include/base/tensor.h index 59c192e4b3..9ad9e46308 100644 --- a/deal.II/base/include/base/tensor.h +++ b/deal.II/base/include/base/tensor.h @@ -230,10 +230,9 @@ class Tensor void clear (); /** - * Determine an estimate for - * the memory consumption (in - * bytes) of this - * object. + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. */ static unsigned int memory_consumption (); diff --git a/deal.II/base/source/data_out_base.cc b/deal.II/base/source/data_out_base.cc index 2c8cf4dc59..f2fff9e6f2 100644 --- a/deal.II/base/source/data_out_base.cc +++ b/deal.II/base/source/data_out_base.cc @@ -2186,7 +2186,7 @@ void DataOutBase::write_gnuplot (const std::vector > &patche { const unsigned int n_subdivisions = patch->n_subdivisions; const unsigned int n = n_subdivisions+1; - // Length of loops in all dimensons + // Length of loops in all dimensions const unsigned int n1 = (dim>0) ? n : 1; const unsigned int n2 = (dim>1) ? n : 1; const unsigned int n3 = (dim>2) ? n : 1; diff --git a/deal.II/base/source/geometry_info.cc b/deal.II/base/source/geometry_info.cc index 2b5ffb638d..646b779dc7 100644 --- a/deal.II/base/source/geometry_info.cc +++ b/deal.II/base/source/geometry_info.cc @@ -17,6 +17,14 @@ DEAL_II_NAMESPACE_OPEN +unsigned int +GeometryInfo<0>::n_children(const RefinementCase<0> &) +{ + return 0; +} + + + template const unsigned int GeometryInfo::max_children_per_cell; template const unsigned int GeometryInfo::faces_per_cell; template const unsigned int GeometryInfo::max_children_per_face; diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 72fd248008..dbd291c40e 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -627,6 +627,319 @@ QWeddle<1>::QWeddle () } +template <> +QGaussLog<1>::QGaussLog(const unsigned int n, + const bool revert) + : + Quadrature<1> (n) +{ + + std::vector p=set_quadrature_points(n); + std::vector w=set_quadrature_weights(n); + + if (revert == true) + revert_points_and_weights(p,w); + + for (unsigned int i=0; isize(); ++i) + { + this->quadrature_points[i] = Point<1>(p[i]); + this->weights[i] = w[i]; + } + +} + +template <> +std::vector +QGaussLog<1>::set_quadrature_points(const unsigned int n) const +{ + + std::vector points(n); + + switch (n) + { + case 1: + points[0] = 0.3333333333333333; + break; + + case 2: + points[0] = 0.1120088061669761; + points[1] = 0.6022769081187381; + break; + + case 3: + points[0] = 0.06389079308732544; + points[1] = 0.3689970637156184; + points[2] = 0.766880303938942; + break; + + case 4: + points[0] = 0.04144848019938324; + points[1] = 0.2452749143206022; + points[2] = 0.5561654535602751; + points[3] = 0.848982394532986; + break; + + case 5: + points[0] = 0.02913447215197205; + points[1] = 0.1739772133208974; + points[2] = 0.4117025202849029; + points[3] = 0.6773141745828183; + points[4] = 0.89477136103101; + break; + + case 6: + points[0] = 0.02163400584411693; + points[1] = 0.1295833911549506; + points[2] = 0.3140204499147661; + points[3] = 0.5386572173517997; + points[4] = 0.7569153373774084; + points[5] = 0.922668851372116; + break; + + + case 7: + points[0] = 0.0167193554082585; + points[1] = 0.100185677915675; + points[2] = 0.2462942462079286; + points[3] = 0.4334634932570557; + points[4] = 0.6323509880476823; + points[5] = 0.81111862674023; + points[6] = 0.940848166743287; + break; + + case 8: + points[0] = 0.01332024416089244; + points[1] = 0.07975042901389491; + points[2] = 0.1978710293261864; + points[3] = 0.354153994351925; + points[4] = 0.5294585752348643; + points[5] = 0.7018145299391673; + points[6] = 0.849379320441094; + points[7] = 0.953326450056343; + break; + + case 9: + points[0] = 0.01086933608417545; + points[1] = 0.06498366633800794; + points[2] = 0.1622293980238825; + points[3] = 0.2937499039716641; + points[4] = 0.4466318819056009; + points[5] = 0.6054816627755208; + points[6] = 0.7541101371585467; + points[7] = 0.877265828834263; + points[8] = 0.96225055941096; + break; + + case 10: + points[0] = 0.00904263096219963; + points[1] = 0.05397126622250072; + points[2] = 0.1353118246392511; + points[3] = 0.2470524162871565; + points[4] = 0.3802125396092744; + points[5] = 0.5237923179723384; + points[6] = 0.6657752055148032; + points[7] = 0.7941904160147613; + points[8] = 0.898161091216429; + points[9] = 0.9688479887196; + break; + + + case 11: + points[0] = 0.007643941174637681; + points[1] = 0.04554182825657903; + points[2] = 0.1145222974551244; + points[3] = 0.2103785812270227; + points[4] = 0.3266955532217897; + points[5] = 0.4554532469286375; + points[6] = 0.5876483563573721; + points[7] = 0.7139638500230458; + points[8] = 0.825453217777127; + points[9] = 0.914193921640008; + points[10] = 0.973860256264123; + break; + + case 12: + points[0] = 0.006548722279080035; + points[1] = 0.03894680956045022; + points[2] = 0.0981502631060046; + points[3] = 0.1811385815906331; + points[4] = 0.2832200676673157; + points[5] = 0.398434435164983; + points[6] = 0.5199526267791299; + points[7] = 0.6405109167754819; + points[8] = 0.7528650118926111; + points[9] = 0.850240024421055; + points[10] = 0.926749682988251; + points[11] = 0.977756129778486; + break; + + default: + Assert(false, ExcNotImplemented()); + break; + } + + return points; +} + + +template <> +std::vector +QGaussLog<1>::set_quadrature_weights(const unsigned int n) const +{ + + std::vector weights(n); + + switch (n) + { + case 1: + weights[0] = -1.0; + break; + case 2: + weights[0] = -0.7185393190303845; + weights[1] = -0.2814606809696154; + break; + + case 3: + weights[0] = -0.5134045522323634; + weights[1] = -0.3919800412014877; + weights[2] = -0.0946154065661483; + break; + + case 4: + weights[0] =-0.3834640681451353; + weights[1] =-0.3868753177747627; + weights[2] =-0.1904351269501432; + weights[3] =-0.03922548712995894; + break; + + case 5: + weights[0] =-0.2978934717828955; + weights[1] =-0.3497762265132236; + weights[2] =-0.234488290044052; + weights[3] =-0.0989304595166356; + weights[4] =-0.01891155214319462; + break; + + case 6: + weights[0] = -0.2387636625785478; + weights[1] = -0.3082865732739458; + weights[2] = -0.2453174265632108; + weights[3] = -0.1420087565664786; + weights[4] = -0.05545462232488041; + weights[5] = -0.01016895869293513; + break; + + + case 7: + weights[0] = -0.1961693894252476; + weights[1] = -0.2703026442472726; + weights[2] = -0.239681873007687; + weights[3] = -0.1657757748104267; + weights[4] = -0.0889432271377365; + weights[5] = -0.03319430435645653; + weights[6] = -0.005932787015162054; + break; + + case 8: + weights[0] = -0.164416604728002; + weights[1] = -0.2375256100233057; + weights[2] = -0.2268419844319134; + weights[3] = -0.1757540790060772; + weights[4] = -0.1129240302467932; + weights[5] = -0.05787221071771947; + weights[6] = -0.02097907374214317; + weights[7] = -0.003686407104036044; + break; + + case 9: + weights[0] = -0.1400684387481339; + weights[1] = -0.2097722052010308; + weights[2] = -0.211427149896601; + weights[3] = -0.1771562339380667; + weights[4] = -0.1277992280331758; + weights[5] = -0.07847890261203835; + weights[6] = -0.0390225049841783; + weights[7] = -0.01386729555074604; + weights[8] = -0.002408041036090773; + break; + + case 10: + weights[0] = -0.12095513195457; + weights[1] = -0.1863635425640733; + weights[2] = -0.1956608732777627; + weights[3] = -0.1735771421828997; + weights[4] = -0.135695672995467; + weights[5] = -0.0936467585378491; + weights[6] = -0.05578772735275126; + weights[7] = -0.02715981089692378; + weights[8] = -0.00951518260454442; + weights[9] = -0.001638157633217673; + break; + + + case 11: + weights[0] = -0.1056522560990997; + weights[1] = -0.1665716806006314; + weights[2] = -0.1805632182877528; + weights[3] = -0.1672787367737502; + weights[4] = -0.1386970574017174; + weights[5] = -0.1038334333650771; + weights[6] = -0.06953669788988512; + weights[7] = -0.04054160079499477; + weights[8] = -0.01943540249522013; + weights[9] = -0.006737429326043388; + weights[10] = -0.001152486965101561; + break; + + case 12: + weights[0] = -0.09319269144393; + weights[1] = -0.1497518275763289; + weights[2] = -0.166557454364573; + weights[3] = -0.1596335594369941; + weights[4] = -0.1384248318647479; + weights[5] = -0.1100165706360573; + weights[6] = -0.07996182177673273; + weights[7] = -0.0524069547809709; + weights[8] = -0.03007108900074863; + weights[9] = -0.01424924540252916; + weights[10] = -0.004899924710875609; + weights[11] = -0.000834029009809656; + break; + + default: + Assert(false, ExcNotImplemented()); + break; + } + + return weights; + +} + + +template <> +void +QGaussLog<1>::revert_points_and_weights(std::vector &p, + std::vector &w) const +{ + + double temp; + const unsigned int sz = this->size(); + + for (unsigned int i=0; i class FullMatrix; template class SparseMatrix; template class Vector; -template class DoFObjectAccessor; -template class DoFObjectAccessor<0,DH>; -template class DoFObjectAccessor<1,DH>; -template class DoFObjectAccessor<2,DH>; -template class DoFObjectAccessor<3,DH>; - - -template class TriaRawIterator; +template class TriaRawIterator; @@ -44,60 +37,56 @@ template class TriaRawIterator; // inlining and thus faster code. +namespace internal +{ + namespace DoFAccessor + { /** * This is a switch class which only declares a @p typedef. It is meant to - * determine which class a DoFObjectAccessor class is to be derived - * from. By default, DoFObjectAccessor derives from - * the @p typedef in the general DoFObjectAccessor_Inheritance - * class, which is TriaObjectAccessor, - * but if celldim==dim, then the specialization DoFObjectAccessor_Inheritance - * is used which declares its local type to be CellAccessor. Therefore, - * the inheritance is automatically chosen to be from CellAccessor if the - * object under consideration has full dimension, i.e. constitutes a cell. + * determine which class a DoFAccessor class is to be derived from. By + * default, DoFAccessor@ derives from the + * typedef in the general + * Inheritance@ class, which is + * TriaAccessor@, but if + * celldim==dim, then the specialization + * Inheritance@ is used which declares + * its local type to be CellAccessor@. Therefore, the + * inheritance is automatically chosen to be from CellAccessor if the object + * under consideration has full dimension, i.e. constitutes a cell. * * @ingroup dofs * @ingroup Accessors * @author Wolfgang Bangerth, 1999 */ -template -class DoFObjectAccessor_Inheritance -{ - public: - /** - * Declaration of the @p typedef. - * See the full documentation for - * more information. - */ - typedef TriaObjectAccessor BaseClass; -}; + template + struct Inheritance + { + /** + * Declaration of the @p typedef. + * See the full documentation for + * more information. + */ + typedef TriaAccessor BaseClass; + }; /** - * This is a switch class which only declares a @p typedef. It is meant to - * determine which class a DoFObjectAccessor class is to be derived - * from. By default, DoFObjectAccessor derives from - * the @p typedef in the general DoFObjectAccessor_Inheritance - * class, which is TriaObjectAccessor, - * but if celldim==dim, then the specialization DoFObjectAccessor_Inheritance - * is used which declares its local type to be CellAccessor. Therefore, - * the inheritance is automatically chosen to be from CellAccessor if the - * object under consideration has full dimension, i.e. constitutes a cell. - * - * @ingroup dofs - * @ingroup Accessors - * @author Wolfgang Bangerth, 1999 + * This is the specialization of the general template used for the case where + * an object has full dimension, i.e. is a cell. See the general template for + * more details. */ -template -class DoFObjectAccessor_Inheritance -{ - public: - /** - * Declaration of the @p typedef. - * See the full documentation for - * more information. - */ - typedef CellAccessor BaseClass; -}; + template + struct Inheritance + { + /** + * Declaration of the @p typedef. + * See the full documentation for + * more information. + */ + typedef CellAccessor BaseClass; + }; + } +} /* -------------------------------------------------------------------------- */ @@ -105,29 +94,26 @@ class DoFObjectAccessor_Inheritance /** - * Define the base class for accessors to the degrees of - * freedom. Accessors are used to, well, access the data that pertains - * to edges, faces, and cells of a triangulation. The concept is + * A class that gives access to the degrees of freedom stored in a DoFHandler + * or hp::DoFHandler object. Accessors are used to, well, access the data that + * pertains to edges, faces, and cells of a triangulation. The concept is * explained in more detail in connection to @ref Iterators. * * This class follows mainly the route laid out by the accessor library - * declared in the triangulation library (TriaAccessor). It enables - * the user to access the degrees of freedom on the lines (there are similar - * versions for the DoFs on quads, etc), where the dimension of the underlying - * triangulation does not really matter (i.e. this accessor works with the - * lines in 1D-, 2D-, etc dimensions). - * - * The first template argument of this class refers to the structural - * dimension of the thing accessed: it is 1 for lines, 2 for quads, or - * 3 for hexes. The second argument denotes the type of DoF handler we - * should work on. It can either be ::DoFHandler or - * hp::DoFHandler. The space dimension of the object we are to - * work on is automatically extracted from the DH template argument. + * declared in the triangulation library (TriaAccessor). It enables the user + * to access the degrees of freedom on lines, quads, or hexes. The first + * template argument of this class determines the dimensionality of the object + * under consideration: 1 for lines, 2 for quads, and 3 for hexes. From the + * second template argument we can deduce the dimensionality of the + * triangulation to which this object belongs as well as the dimensionality of + * the space in which it is embedded. The second argument also denotes the + * type of DoF handler we should work on. It can either be + * ::DoFHandler or hp::DoFHandler. * * Depending on whether the structural dimension of the object - * accessed equals the space dimension on which the DoF handler object + * accessed equals the dimension on which the DoF handler object * operates, this class is derived from CellAccessor or - * TriaObjectAccessor. This means that, for example accessors to quads + * TriaAccessor. This means that, for example accessors to quads * in 2d have access to all the mesh aspects of cells, whereas * accessors to quads in 3d can only access things that make sense for * faces. @@ -150,28 +136,34 @@ class DoFObjectAccessor_Inheritance * rather than lower-dimensional objects. In that case, inheritance is * from CellAccessor, to provide access to all the cell specific * information afforded by that class. Otherwise, i.e. for - * lower-dimensional objects, inheritance is from TriaObjectAccessor. - * - * There are classes DoFObjectAccessor<1>, DoFObjectAccessor<2>, and - * DoFObjectAccessor<3> derived from the present class that provide - * things that are specific to objects of a particular - * dimensionality. For example, DoFObjectAccessor<2> allows access to - * the DoF iterators pointing to the lines that bound that - * quadrilateral, which is obviously something that a line can't - * do. Consequently, these few operations are declared in these - * derived classes. + * lower-dimensional objects, inheritance is from TriaAccessor. * - * In addition, there is a DoFCellAccessor class that provides the + * There is a DoFCellAccessor class that provides the * equivalent to the CellAccessor class. * * @ingroup dofs * @ingroup Accessors - * @author Wolfgang Bangerth, 1998, 2006 + * @author Wolfgang Bangerth, 1998, 2006, 2008 */ template -class DoFAccessor : public DoFObjectAccessor_Inheritance::BaseClass +class DoFAccessor : public internal::DoFAccessor::Inheritance::BaseClass { public: + + /** + * A static variable that allows users of + * this class to discover the value of + * the second template argument. + */ + static const unsigned int dimension=DH::dimension; + + /** + * A static variable that allows users of + * this class to discover the value of + * the third template argument. + */ + static const unsigned int space_dimension=DH::space_dimension; + /** * Declare a typedef to the base * class to make accessing some @@ -179,7 +171,7 @@ class DoFAccessor : public DoFObjectAccessor_Inheritance::BaseClass + typename internal::DoFAccessor::Inheritance::BaseClass BaseClass; /** @@ -187,6 +179,13 @@ class DoFAccessor : public DoFObjectAccessor_Inheritance *tria, + DoFAccessor (const Triangulation *tria, const int level, const int index, const DH *local_data); /** - * Reset the DoF handler pointer. + * Conversion constructor. This + * constructor exists to make certain + * constructs simpler to write in + * dimension independent code. For + * example, it allows assigning a face + * iterator to a line iterator, an + * operation that is useful in 2d but + * doesn't make any sense in 3d. The + * constructor here exists for the + * purpose of making the code conform to + * C++ but it will unconditionally abort; + * in other words, assigning a face + * iterator to a line iterator is better + * put into an if-statement that checks + * that the dimension is two, and assign + * to a quad iterator in 3d (an operator + * that, without this constructor would + * be illegal if we happen to compile for + * 2d). + */ + template + DoFAccessor (const InvalidAccessor &); + + /** + * Another conversion operator + * between objects that don't + * make sense, just like the + * previous one. */ - void set_dof_handler (DH *dh); + template + DoFAccessor (const DoFAccessor &); + /** + * @} + */ + /** * Return a handle on the * DoFHandler object which we @@ -227,17 +258,71 @@ class DoFAccessor : public DoFObjectAccessor_Inheritance &a); + /** + * @name Accessing sub-objects + */ + /** + * @{ + */ + /** * Return an iterator pointing to * the the @p c-th child. */ - TriaIterator > + TriaIterator > child (const unsigned int c) const; /** - * Global DoF index of the i - * degree associated with the @p vertexth - * vertex of the present cell. + * Pointer to the @p ith line + * bounding this object. + */ + typename internal::DoFHandler::Iterators::line_iterator + line (const unsigned int i) const; + + /** + * Pointer to the @p ith quad + * bounding this object. + */ + typename internal::DoFHandler::Iterators::quad_iterator + quad (const unsigned int i) const; + + /** + * @} + */ + + /** + * @name Accessing the DoF indices of this object + */ + /** + * @{ + */ + + /** + * Return the indices of the dofs of this + * object in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, etc, + * dofs on line 0, dofs on line 1, etc, + * dofs on quad 0, etc. + * + * The vector has to have the + * right size before being passed + * to this function. + * + * This function is most often + * used on active objects (edges, + * faces, cells). It can be used + * on non-active objects as well + * (i.e. objects that have + * children), but only if the + * finite element under + * consideration has degrees of + * freedom exclusively on + * vertices. Otherwise, the + * function doesn't make much + * sense, since for example + * inactive edges do not have + * degrees of freedom associated + * with them at all. * * The last argument denotes the * finite element index. For the @@ -269,15 +354,27 @@ class DoFAccessor : public DoFObjectAccessor_Inheritancecell-@>active_fe_index. Consequently, + * the derived DoFCellAccessor + * class has an overloaded + * version of this function that + * calls the present function + * with + * cell-@>active_fe_index + * as last argument. */ - unsigned int vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const unsigned int fe_index = DH::default_fe_index) const; + void get_dof_indices (std::vector &dof_indices, + const unsigned int fe_index = DH::default_fe_index) const; /** - * Set the global index of the i - * degree on the @p vertex-th vertex of - * the present cell to @p index. + * Global DoF index of the i + * degree associated with the @p vertexth + * vertex of the present cell. * * The last argument denotes the * finite element index. For the @@ -310,10 +407,9 @@ class DoFAccessor : public DoFObjectAccessor_Inheritanceith degree @@ -353,46 +449,17 @@ class DoFAccessor : public DoFObjectAccessor_Inheritanceith degree of freedom - * of this object to @p index. - * - * The last argument denotes the - * finite element index. For the - * standard ::DoFHandler class, - * this value must be equal to - * its default value since that - * class only supports the same - * finite element on all cells - * anyway. - * - * However, for hp objects - * (i.e. the hp::DoFHandler - * class), different finite - * element objects may be used on - * different cells. On faces - * between two cells, as well as - * vertices, there may therefore - * be two sets of degrees of - * freedom, one for each of the - * finite elements used on the - * adjacent cells. In order to - * specify which set of degrees - * of freedom to work on, the - * last argument is used to - * disambiguate. Finally, if this - * function is called for a cell - * object, there can only be a - * single set of degrees of - * freedom, and fe_index has to - * match the result of - * active_fe_index(). + /** + * @} */ - void set_dof_index (const unsigned int i, - const unsigned int index, - const unsigned int fe_index = DH::default_fe_index) const; + /** + * @name Accessing the finite element associated with this object + */ + /** + * @{ + */ + /** * Return the number of finite * elements that are active on a @@ -461,8 +528,12 @@ class DoFAccessor : public DoFObjectAccessor_Inheritancefe_index_is_active(fe_index) * must return true. */ - const FiniteElement & + const FiniteElement & get_fe (const unsigned int fe_index) const; + + /** + * @} + */ /** * Exceptions for child classes @@ -505,261 +576,31 @@ class DoFAccessor : public DoFObjectAccessor_Inheritance friend class TriaRawIterator; -}; - - -/* -------------------------------------------------------------------------- */ - - -/** - * Closure class. Unused, but necessary for some template constructs. - * @ingroup dofs - * @ingroup Accessors - */ -template -class DoFObjectAccessor<0, DH> : public DoFAccessor<0, DH> -{ - public: - /** - * Extract dimension from DH. - */ - static const unsigned int dim = DH::dimension; - - typedef void AccessorData; - - /** - * Constructor. Should never be called - * and thus throws an error. - */ - DoFObjectAccessor (Triangulation *, - const int, - const int, - const AccessorData *) - { - Assert (false, ExcInternalError()); - } -}; - - - -/** - * Grant access to those aspects of degrees of freedom located on - * quads that are dimension specific. See the DoFAccessor class for - * more information. - * - * @ingroup dofs - * @ingroup Accessors - * @author Wolfgang Bangerth, 1998, 2006 - */ -template -class DoFObjectAccessor<1, DH> : public DoFAccessor<1,DH> -{ - public: - /** - * Extract dimension from DH. - */ - static const unsigned int dim = DH::dimension; - - /** - * Declare the data type that this accessor - * class expects to get passed from the - * iterator classes. - */ - typedef DH AccessorData; - - /** - * Declare base class as a local typedef - * for simpler access. - */ - typedef DoFAccessor<1,DH> BaseClass; - - /** - * Default constructor, unused thus - * not implemented. - */ - DoFObjectAccessor (); - - /** - * Constructor. The @p local_data - * argument is assumed to be a pointer - * to a DoFHandler object. - */ - DoFObjectAccessor (const Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data); - - /** - * Return the indices of the dofs of this - * object in the standard ordering: dofs - * on vertex 0, dofs on vertex 1, etc, - * dofs on line 0, dofs on line 1, etc, - * dofs on quad 0, etc. - * - * The vector has to have the - * right size before being passed - * to this function. - * - * This function is most often - * used on active objects (edges, - * faces, cells). It can be used - * on non-active objects as well - * (i.e. objects that have - * children), but only if the - * finite element under - * consideration has degrees of - * freedom exclusively on - * vertices. Otherwise, the - * function doesn't make much - * sense, since for example - * inactive edges do not have - * degrees of freedom associated - * with them at all. - * - * The last argument denotes the - * finite element index. For the - * standard ::DoFHandler class, - * this value must be equal to - * its default value since that - * class only supports the same - * finite element on all cells - * anyway. - * - * However, for hp objects - * (i.e. the hp::DoFHandler - * class), different finite - * element objects may be used on - * different cells. On faces - * between two cells, as well as - * vertices, there may therefore - * be two sets of degrees of - * freedom, one for each of the - * finite elements used on the - * adjacent cells. In order to - * specify which set of degrees - * of freedom to work on, the - * last argument is used to - * disambiguate. Finally, if this - * function is called for a cell - * object, there can only be a - * single set of degrees of - * freedom, and fe_index has to - * match the result of - * active_fe_index(). - * - * For cells, there is only a - * single possible finite element - * index (namely the one for that - * cell, returned by - * cell-@>active_fe_index. Consequently, - * the derived DoFCellAccessor - * class has an overloaded - * version of this function that - * calls the present function - * with - * cell-@>active_fe_index - * as last argument. - */ - void get_dof_indices (std::vector &dof_indices, - const unsigned int fe_index = DH::default_fe_index) const; -}; - - -/** - * Grant access to those aspects of degrees of freedom located on - * quads that are dimension specific. See the DoFAccessor class for - * more information. - * - * @ingroup dofs - * @ingroup Accessors - * @author Wolfgang Bangerth, 1998, 2006 - */ -template -class DoFObjectAccessor<2, DH> : public DoFAccessor<2,DH> -{ - public: - /** - * Extract dimension from DH. - */ - static const unsigned int dim = DH::dimension; - - /** - * Declare the data type that this accessor - * class expects to get passed from the - * iterator classes. - */ - typedef DH AccessorData; - - /** - * Declare base class as a local typedef - * for simpler access. - */ - typedef DoFAccessor<2,DH> BaseClass; - - /** - * Default constructor, unused thus - * not implemented. - */ - DoFObjectAccessor (); - - /** - * Constructor. The @p local_data - * argument is assumed to be a pointer - * to a DoFHandler object. - */ - DoFObjectAccessor (const Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data); - - /** - * Return the indices of the dofs of this - * object in the standard ordering: dofs - * on vertex 0, dofs on vertex 1, etc, - * dofs on line 0, dofs on line 1, etc, - * dofs on quad 0, etc. - * - * The vector has to have the - * right size before being passed - * to this function. - * - * This function is most often - * used on active objects (edges, - * faces, cells). It can be used - * on non-active objects as well - * (i.e. objects that have - * children), but only if the - * finite element under - * consideration has degrees of - * freedom exclusively on - * vertices. Otherwise, the - * function doesn't make much - * sense, since for example - * inactive edges do not have - * degrees of freedom associated - * with them at all. + /** + * Store the address of the DoFHandler object + * to be accessed. + */ + DH *dof_handler; + + /** + * Compare for equality. + */ + bool operator == (const DoFAccessor &) const; + + /** + * Compare for inequality. + */ + bool operator != (const DoFAccessor &) const; + + /** + * Reset the DoF handler pointer. + */ + void set_dof_handler (DH *dh); + + /** + * Set the index of the + * ith degree of freedom + * of this object to @p index. * * The last argument denotes the * finite element index. For the @@ -791,105 +632,15 @@ class DoFObjectAccessor<2, DH> : public DoFAccessor<2,DH> * freedom, and fe_index has to * match the result of * active_fe_index(). - * - * For cells, there is only a - * single possible finite element - * index (namely the one for that - * cell, returned by - * cell-@>active_fe_index. Consequently, - * the derived DoFCellAccessor - * class has an overloaded - * version of this function that - * calls the present function - * with - * cell-@>active_fe_index - * as last argument. - */ - void get_dof_indices (std::vector &dof_indices, - const unsigned int fe_index = DH::default_fe_index) const; - - /** - * Return a pointer to the @p ith line - * bounding this @p Quad. - */ - TriaIterator > - line (const unsigned int i) const; -}; - - -/** - * Grant access to those aspects of degrees of freedom located on - * hexes that are dimension specific. See the DoFAccessor class for - * more information. - * - * @ingroup dofs - * @ingroup Accessors - * @author Wolfgang Bangerth, 1998, 2006 - */ -template -class DoFObjectAccessor<3, DH> : public DoFAccessor<3,DH> -{ - public: - /** - * Extract dimension from DH. - */ - static const unsigned int dim = DH::dimension; - - /** - * Declare the data type that this accessor - * class expects to get passed from the - * iterator classes. */ - typedef DH AccessorData; + void set_dof_index (const unsigned int i, + const unsigned int index, + const unsigned int fe_index = DH::default_fe_index) const; /** - * Declare base class as a local typedef - * for simpler access. - */ - typedef DoFAccessor<1,DH> BaseClass; - - /** - * Default constructor, unused thus - * not implemented. - */ - DoFObjectAccessor (); - - /** - * Constructor. The @p local_data - * argument is assumed to be a pointer - * to a DoFHandler object. - */ - DoFObjectAccessor (const Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data); - - /** - * Return the indices of the dofs of this - * object in the standard ordering: dofs - * on vertex 0, dofs on vertex 1, etc, - * dofs on line 0, dofs on line 1, etc, - * dofs on quad 0, etc. - * - * The vector has to have the - * right size before being passed - * to this function. - * - * This function is most often - * used on active objects (edges, - * faces, cells). It can be used - * on non-active objects as well - * (i.e. objects that have - * children), but only if the - * finite element under - * consideration has degrees of - * freedom exclusively on - * vertices. Otherwise, the - * function doesn't make much - * sense, since for example - * inactive edges do not have - * degrees of freedom associated - * with them at all. + * Set the global index of the i + * degree on the @p vertex-th vertex of + * the present cell to @p index. * * The last argument denotes the * finite element index. For the @@ -921,52 +672,37 @@ class DoFObjectAccessor<3, DH> : public DoFAccessor<3,DH> * freedom, and fe_index has to * match the result of * active_fe_index(). - * - * For cells, there is only a - * single possible finite element - * index (namely the one for that - * cell, returned by - * cell-@>active_fe_index. Consequently, - * the derived DoFCellAccessor - * class has an overloaded - * version of this function that - * calls the present function - * with - * cell-@>active_fe_index - * as last argument. */ - void get_dof_indices (std::vector &dof_indices, - const unsigned int fe_index = DH::default_fe_index) const; + void set_vertex_dof_index (const unsigned int vertex, + const unsigned int i, + const unsigned int index, + const unsigned int fe_index = DH::default_fe_index) const; + + /** + * Iterator classes need to be friends + * because they need to access operator== + * and operator!=. + */ + template friend class TriaRawIterator; - /** - * Return a pointer to the @p ith line - * bounding this @p Hex. - */ - TriaIterator > - line (const unsigned int i) const; - /** - * Return a pointer to the @p ith quad - * bounding this @p Hex. + /** + * Make the DoFHandler class a friend so + * that it can call the set_xxx() + * functions. */ - TriaIterator > - quad (const unsigned int i) const; + template friend class DoFHandler; + template friend class hp::DoFHandler; }; +/* -------------------------------------------------------------------------- */ + /** - * Grant access to the degrees of freedom on a cell. In fact, since all - * access to the degrees of freedom has been enabled by the classes - * DoFObjectAccessor<1, 1> and DoFObjectAccessor<2, 2> for the space dimension - * one and two, respectively, this class only collects the pieces - * together by deriving from the appropriate DoF*Accessor and the - * right CellAccessor and finally adding two functions which give - * access to the neighbors and children as DoFCellAccessor objects - * rather than CellAccessor objects (the latter function was inherited - * from the CellAccessor class). + * Grant access to the degrees of freedom on a cell. * - * Note that since for the class we derive from, i.e. DoFObjectAccessor, + * Note that since for the class we derive from, i.e. DoFAccessor, * the two template parameters are equal, the base class is actually derived from * CellAccessor, which makes the functions of this class available to the * DoFCellAccessor class as well. @@ -976,13 +712,18 @@ class DoFObjectAccessor<3, DH> : public DoFAccessor<3,DH> * @author Wolfgang Bangerth, 1998 */ template -class DoFCellAccessor : public DoFObjectAccessor +class DoFCellAccessor : public DoFAccessor { public: /** * Extract dimension from DH. */ static const unsigned int dim = DH::dimension; + + /** + * Extract space dimension from DH. + */ + static const unsigned int spacedim = DH::space_dimension; /** * Declare the data type that @@ -990,7 +731,7 @@ class DoFCellAccessor : public DoFObjectAccessor * get passed from the iterator * classes. */ - typedef typename DoFObjectAccessor::AccessorData AccessorData; + typedef typename DoFAccessor::AccessorData AccessorData; /** * Declare a typedef to the base @@ -998,16 +739,67 @@ class DoFCellAccessor : public DoFObjectAccessor * of the exception classes * simpler. */ - typedef DoFObjectAccessor BaseClass; + typedef DoFAccessor BaseClass; + /** + * @name Constructors + */ + /** + * @{ + */ + /** * Constructor */ - DoFCellAccessor (const Triangulation *tria, + DoFCellAccessor (const Triangulation *tria, const int level, const int index, const AccessorData *local_data); + /** + * Conversion constructor. This + * constructor exists to make certain + * constructs simpler to write in + * dimension independent code. For + * example, it allows assigning a face + * iterator to a line iterator, an + * operation that is useful in 2d but + * doesn't make any sense in 3d. The + * constructor here exists for the + * purpose of making the code conform to + * C++ but it will unconditionally abort; + * in other words, assigning a face + * iterator to a line iterator is better + * put into an if-statement that checks + * that the dimension is two, and assign + * to a quad iterator in 3d (an operator + * that, without this constructor would + * be illegal if we happen to compile for + * 2d). + */ + template + DoFCellAccessor (const InvalidAccessor &); + + /** + * Another conversion operator + * between objects that don't + * make sense, just like the + * previous one. + */ + template + DoFCellAccessor (const DoFAccessor &); + + /** + * @} + */ + + /** + * @name Accessing sub-objects and neighbors + */ + /** + * @{ + */ + /** * Return the @p ith neighbor as * a DoF cell iterator. This @@ -1017,7 +809,7 @@ class DoFCellAccessor : public DoFObjectAccessor * without access to the DoF * data. */ - TriaIterator > + typename internal::DoFHandler::Iterators::cell_iterator neighbor (const unsigned int) const; /** @@ -1029,7 +821,7 @@ class DoFCellAccessor : public DoFObjectAccessor * without access to the DoF * data. */ - TriaIterator > + typename internal::DoFHandler::Iterators::cell_iterator child (const unsigned int) const; /** @@ -1037,10 +829,10 @@ class DoFCellAccessor : public DoFObjectAccessor * of this cell. * * This function is not implemented in - * 1D, and maps to DoFObjectAccessor<2, + * 1D, and maps to DoFAccessor<2, * dim>::line in 2D. */ - TriaIterator > + typename internal::DoFHandler::Iterators::face_iterator face (const unsigned int i) const; /** @@ -1054,10 +846,21 @@ class DoFCellAccessor : public DoFObjectAccessor * access to the triangulation * data). */ - TriaIterator > + typename internal::DoFHandler::Iterators::cell_iterator neighbor_child_on_subface (const unsigned int face_no, const unsigned int subface_no) const; + /** + * @} + */ + + /** + * @name Extracting values from global vectors + */ + /** + * @{ + */ + /** * Return the values of the given vector * restricted to the dofs of this @@ -1257,6 +1060,54 @@ class DoFCellAccessor : public DoFObjectAccessor void set_dof_values_by_interpolation (const Vector &local_values, OutputVector &values) const; + /** + * Distribute a local (cell + * based) vector to a global one + * by mapping the local numbering + * of the degrees of freedom to + * the global one and entering + * the local values into the + * global vector. + * + * The elements are + * added up to the + * elements in the global vector, + * rather than just set, since + * this is usually what one + * wants. + */ + template + void + distribute_local_to_global (const Vector &local_source, + OutputVector &global_destination) const; + + /** + * This function does much the + * same as the + * distribute_local_to_global(Vector,Vector) + * function, but operates on + * matrices instead of + * vectors. If the matrix type is + * a sparse matrix then it is + * supposed to have non-zero + * entry slots where required. + */ + template + void + distribute_local_to_global (const FullMatrix &local_source, + OutputMatrix &global_destination) const; + + /** + * @} + */ + + /** + * @name Accessing the DoF indices of this object + */ + /** + * @{ + */ + /** * Return the indices of the dofs of this * quad in the standard ordering: dofs @@ -1304,6 +1155,17 @@ class DoFCellAccessor : public DoFObjectAccessor */ void get_dof_indices (std::vector &dof_indices) const; + /** + * @} + */ + + /** + * @name Accessing the finite element associated with this object + */ + /** + * @{ + */ + /** * Return the finite element that * is used on the cell pointed to @@ -1315,7 +1177,7 @@ class DoFCellAccessor : public DoFObjectAccessor * for hp DoF handlers, this may * change from cell to cell. */ - const FiniteElement & + const FiniteElement & get_fe () const; /** @@ -1330,45 +1192,21 @@ class DoFCellAccessor : public DoFObjectAccessor * this cell. */ void set_active_fe_index (const unsigned int i); - /** - * Distribute a local (cell - * based) vector to a global one - * by mapping the local numbering - * of the degrees of freedom to - * the global one and entering - * the local values into the - * global vector. - * - * The elements are - * added up to the - * elements in the global vector, - * rather than just set, since - * this is usually what one - * wants. + * @} */ - template - void - distribute_local_to_global (const Vector &local_source, - OutputVector &global_destination) const; - + + private: /** - * This function does much the - * same as the - * distribute_local_to_global(Vector,Vector) - * function, but operates on - * matrices instead of - * vectors. If the matrix type is - * a sparse matrix then it is - * supposed to have non-zero - * entry slots where required. + * Forward declaration of a class + * into which we put some + * parts of the implementation. + * + * See the .cc file for more + * information. */ - template - void - distribute_local_to_global (const FullMatrix &local_source, - OutputMatrix &global_destination) const; + struct Implementation; - private: /** * Update the cache in which we * store the dof indices of this @@ -1382,74 +1220,10 @@ class DoFCellAccessor : public DoFObjectAccessor * update_cell_dof_indices_cache() * function */ - template friend class DoFHandler; + template friend class DoFHandler; }; -/* -------------- declaration of explicit specializations ------------- */ - -template <> -TriaIterator<1, DoFObjectAccessor<0,DoFHandler<1> > > -DoFCellAccessor >::face (const unsigned int) const; -template <> -TriaIterator<2, DoFObjectAccessor<1,DoFHandler<2> > > -DoFCellAccessor >::face (const unsigned int i) const; -template <> -TriaIterator<3, DoFObjectAccessor<2,DoFHandler<3> > > -DoFCellAccessor >::face (const unsigned int i) const; - - -template <> -const FiniteElement<1> & -DoFCellAccessor >::get_fe () const; -template <> -const FiniteElement<2> & -DoFCellAccessor >::get_fe () const; -template <> -const FiniteElement<3> & -DoFCellAccessor >::get_fe () const; - -template <> -unsigned int -DoFCellAccessor >::active_fe_index () const; -template <> -unsigned int -DoFCellAccessor >::active_fe_index () const; -template <> -unsigned int -DoFCellAccessor >::active_fe_index () const; - -template <> -void -DoFCellAccessor >::set_active_fe_index (const unsigned int i); -template <> -void -DoFCellAccessor >::set_active_fe_index (const unsigned int i); -template <> -void -DoFCellAccessor >::set_active_fe_index (const unsigned int i); - -template <> -void -DoFCellAccessor >::update_cell_dof_indices_cache () const; -template <> -void -DoFCellAccessor >::update_cell_dof_indices_cache () const; -template <> -void -DoFCellAccessor >::update_cell_dof_indices_cache () const; - -template <> -void -DoFCellAccessor >::update_cell_dof_indices_cache () const; -template <> -void -DoFCellAccessor >::update_cell_dof_indices_cache () const; -template <> -void -DoFCellAccessor >::update_cell_dof_indices_cache () const; - - DEAL_II_NAMESPACE_CLOSE // include more templates diff --git a/deal.II/deal.II/include/dofs/dof_accessor.templates.h b/deal.II/deal.II/include/dofs/dof_accessor.templates.h index b78b44a424..2407a84df1 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.templates.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -31,12 +31,8 @@ DEAL_II_NAMESPACE_OPEN template -DoFAccessor::DoFAccessor () - : - DoFObjectAccessor_Inheritance::BaseClass (0, - numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int), - dof_handler(0) +DoFAccessor:: +DoFAccessor () { Assert (false, ExcInvalidObject()); } @@ -45,19 +41,40 @@ DoFAccessor::DoFAccessor () template inline -DoFAccessor::DoFAccessor (const Triangulation *tria, - const int level, - const int index, - const DH *dof_handler) +DoFAccessor:: +DoFAccessor (const Triangulation *tria, + const int level, + const int index, + const DH *dof_handler) : - DoFObjectAccessor_Inheritance::BaseClass (tria, - level, - index), + internal::DoFAccessor::Inheritance::BaseClass (tria, + level, + index), dof_handler(const_cast(dof_handler)) {} +template +template +DoFAccessor:: +DoFAccessor (const InvalidAccessor &) +{ + Assert (false, ExcInvalidObject()); +} + + + +template +template +DoFAccessor:: +DoFAccessor (const DoFAccessor &) +{ + Assert (false, ExcInvalidObject()); +} + + + template inline void @@ -124,7 +141,7 @@ DoFAccessor::operator != (const DoFAccessor &a) cons template inline -TriaIterator > +TriaIterator > DoFAccessor::child (const unsigned int i) const { Assert (static_cast(this->level()) < this->dof_handler->levels.size(), @@ -136,11 +153,11 @@ DoFAccessor::child (const unsigned int i) const else next_level = 0; - TriaIterator > q (this->tria, - next_level, - this->child_index (i), - this->dof_handler); - + TriaIterator > q (this->tria, + next_level, + this->child_index (i), + this->dof_handler); + // make sure that we either created // a past-the-end iterator or one // pointing to a used cell @@ -259,200 +276,200 @@ DoFAccessor::fe_index_is_active (const unsigned int fe_index) const -template -inline -const FiniteElement & -DoFAccessor::get_fe (const unsigned int fe_index) const +namespace internal { - Assert (fe_index_is_active (fe_index) == true, - ExcMessage ("This function can only be called for active fe indices")); - - return this->dof_handler->get_fe()[fe_index]; -} - - - -template <> -inline -const FiniteElement<1> & -DoFAccessor<1,dealii::DoFHandler<1> >::get_fe (const unsigned int fe_index) const -{ - Assert (fe_index_is_active (fe_index) == true, - ExcMessage ("This function can only be called for active fe indices")); - return this->dof_handler->get_fe(); -} - - - -template <> -inline -const FiniteElement<2> & -DoFAccessor<1,dealii::DoFHandler<2> >::get_fe (const unsigned int fe_index) const -{ - Assert (fe_index_is_active (fe_index) == true, - ExcMessage ("This function can only be called for active fe indices")); - return this->dof_handler->get_fe(); -} - + namespace DoFAccessor + { + template + inline + const FiniteElement & + get_fe (const FiniteElement &fe, + const unsigned int) + { + return fe; + } -template <> -inline -const FiniteElement<3> & -DoFAccessor<1,dealii::DoFHandler<3> >::get_fe (const unsigned int fe_index) const -{ - Assert (fe_index_is_active (fe_index) == true, - ExcMessage ("This function can only be called for active fe indices")); - return this->dof_handler->get_fe(); -} - - - -template <> -inline -const FiniteElement<2> & -DoFAccessor<2,dealii::DoFHandler<2> >::get_fe (const unsigned int fe_index) const -{ - Assert (fe_index_is_active (fe_index) == true, - ExcMessage ("This function can only be called for active fe indices")); - return this->dof_handler->get_fe(); + + template + inline + const FiniteElement & + get_fe (const dealii::hp::FECollection &fe, + const unsigned int index) + { + return fe[index]; + } + } } - -template <> +template inline -const FiniteElement<3> & -DoFAccessor<2,dealii::DoFHandler<3> >::get_fe (const unsigned int fe_index) const +const FiniteElement & +DoFAccessor::get_fe (const unsigned int fe_index) const { Assert (fe_index_is_active (fe_index) == true, - ExcMessage ("This function can only be called for active fe indices")); - return this->dof_handler->get_fe(); + ExcMessage ("This function can only be called for active fe indices")); + + return internal::DoFAccessor::get_fe (this->dof_handler->get_fe(), fe_index); } -template <> -inline -const FiniteElement<3> & -DoFAccessor<3,dealii::DoFHandler<3> >::get_fe (const unsigned int fe_index) const +namespace internal { - Assert (fe_index_is_active (fe_index) == true, - ExcMessage ("This function can only be called for active fe indices")); - return this->dof_handler->get_fe(); + namespace DoFAccessor + { + template + void get_dof_indices (const dealii::DoFAccessor<1,DH> &accessor, + std::vector &dof_indices, + const unsigned int fe_index) + { + const unsigned int dofs_per_vertex = accessor.get_fe(fe_index).dofs_per_vertex, + dofs_per_line = accessor.get_fe(fe_index).dofs_per_line; + std::vector::iterator next = dof_indices.begin(); + for (unsigned int vertex=0; vertex<2; ++vertex) + for (unsigned int d=0; d + void get_dof_indices (const dealii::DoFAccessor<2,DH> &accessor, + std::vector &dof_indices, + const unsigned int fe_index) + { + const unsigned int dofs_per_vertex = accessor.get_fe(fe_index).dofs_per_vertex, + dofs_per_line = accessor.get_fe(fe_index).dofs_per_line, + dofs_per_quad = accessor.get_fe(fe_index).dofs_per_quad; + std::vector::iterator next = dof_indices.begin(); + for (unsigned int vertex=0; vertex<4; ++vertex) + for (unsigned int d=0; ddof_index(accessor.get_fe(fe_index). + adjust_line_dof_index_for_line_orientation(d, + accessor.line_orientation(line)), + fe_index); + for (unsigned int d=0; d + void get_dof_indices (const dealii::DoFAccessor<3,DH> &accessor, + std::vector &dof_indices, + const unsigned int fe_index) + { + const unsigned int dofs_per_vertex = accessor.get_fe(fe_index).dofs_per_vertex, + dofs_per_line = accessor.get_fe(fe_index).dofs_per_line, + dofs_per_quad = accessor.get_fe(fe_index).dofs_per_quad, + dofs_per_hex = accessor.get_fe(fe_index).dofs_per_hex; + std::vector::iterator next = dof_indices.begin(); + for (unsigned int vertex=0; vertex<8; ++vertex) + for (unsigned int d=0; ddof_index(accessor.get_fe(fe_index). + adjust_line_dof_index_for_line_orientation(d, + accessor.line_orientation(line)),fe_index); + // now copy dof numbers from the face. for + // faces with the wrong orientation, we + // have already made sure that we're ok by + // picking the correct lines and vertices + // (this happens automatically in the + // line() and vertex() functions). however, + // if the face is in wrong orientation, we + // look at it in flipped orientation and we + // will have to adjust the shape function + // indices that we see to correspond to the + // correct (cell-local) ordering. The same + // applies, if the face_rotation or + // face_orientation is non-standard + for (unsigned int quad=0; quad<6; ++quad) + for (unsigned int d=0; ddof_index(accessor.get_fe(fe_index). + adjust_quad_dof_index_for_face_orientation(d, + accessor.face_orientation(quad), + accessor.face_flip(quad), + accessor.face_rotation(quad)), + fe_index); + for (unsigned int d=0; d -----------------------*/ - -template -inline -DoFObjectAccessor<1,DH>:: -DoFObjectAccessor (const Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data) - : - DoFAccessor<1,DH> (tria, level, index, local_data) -{} - - - - -template +template inline void -DoFObjectAccessor<1,DH>::get_dof_indices (std::vector &dof_indices, - const unsigned int fe_index) const +DoFAccessor::get_dof_indices (std::vector &dof_indices, + const unsigned int fe_index) const { Assert (static_cast(this->level()) < this->dof_handler->levels.size(), ExcMessage ("DoFHandler not initialized")); - Assert (this->dof_handler != 0, typename BaseClass::ExcInvalidObject()); - Assert (&this->dof_handler->get_fe() != 0, typename BaseClass::ExcInvalidObject()); - Assert (dof_indices.size() == (2*this->dof_handler->get_fe()[fe_index].dofs_per_vertex + - this->dof_handler->get_fe()[fe_index].dofs_per_line), - typename BaseClass::ExcVectorDoesNotMatch()); - - // this function really only makes - // sense if either a) there are - // degrees of freedom defined on - // the present object, or b) the - // object is non-active objects but - // all degrees of freedom are - // located on vertices, since - // otherwise there are degrees of - // freedom on sub-objects which are - // not allocated for this - // non-active thing - Assert (this->fe_index_is_active (fe_index) - || - (this->dof_handler->get_fe()[fe_index].dofs_per_cell == - 2*this->dof_handler->get_fe()[fe_index].dofs_per_vertex), - ExcInternalError()); - - const unsigned int dofs_per_vertex = this->dof_handler->get_fe()[fe_index].dofs_per_vertex, - dofs_per_line = this->dof_handler->get_fe()[fe_index].dofs_per_line; - std::vector::iterator next = dof_indices.begin(); - for (unsigned int vertex=0; vertex<2; ++vertex) - for (unsigned int d=0; dvertex_dof_index(vertex,d,fe_index); - for (unsigned int d=0; ddof_index(d,fe_index); -} - - - - - -/*------------------------- Functions: DoFObjectAccessor<2,dim> -----------------------*/ - -template -DoFObjectAccessor<2,DH>:: -DoFObjectAccessor (const Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data) - : - DoFAccessor<2,DH> (tria, level, index, local_data) -{} - - - -template -inline -TriaIterator > -DoFObjectAccessor<2,DH>::line (const unsigned int i) const -{ - Assert (i<4, ExcIndexRange (i, 0, 4)); + Assert (this->dof_handler != 0, ExcInvalidObject()); + Assert (&this->dof_handler->get_fe() != 0, ExcInvalidObject()); Assert (static_cast(this->level()) < this->dof_handler->levels.size(), ExcMessage ("DoFHandler not initialized")); - return TriaIterator > - ( - this->tria, - 0, - this->line_index (i), - this->dof_handler - ); -} - - - -template -inline -void -DoFObjectAccessor<2,DH>::get_dof_indices (std::vector &dof_indices, - const unsigned int fe_index) const -{ - Assert (this->dof_handler != 0, - typename BaseClass::ExcInvalidObject()); - Assert (&this->dof_handler->get_fe() != 0, - typename BaseClass::ExcInvalidObject()); - Assert (dof_indices.size() == (4*this->dof_handler->get_fe()[fe_index].dofs_per_vertex + - 4*this->dof_handler->get_fe()[fe_index].dofs_per_line + - this->dof_handler->get_fe()[fe_index].dofs_per_quad), - typename BaseClass::ExcVectorDoesNotMatch()); + switch (structdim) + { + case 1: + Assert (dof_indices.size() == + (2*this->dof_handler->get_fe()[fe_index].dofs_per_vertex + + this->dof_handler->get_fe()[fe_index].dofs_per_line), + ExcVectorDoesNotMatch()); + break; + case 2: + Assert (dof_indices.size() == + (4*this->dof_handler->get_fe()[fe_index].dofs_per_vertex + + 4*this->dof_handler->get_fe()[fe_index].dofs_per_line + + this->dof_handler->get_fe()[fe_index].dofs_per_quad), + ExcVectorDoesNotMatch()); + break; + case 3: + Assert (dof_indices.size() == + (8*this->dof_handler->get_fe()[fe_index].dofs_per_vertex + + 12*this->dof_handler->get_fe()[fe_index].dofs_per_line + + 6*this->dof_handler->get_fe()[fe_index].dofs_per_quad + + this->dof_handler->get_fe()[fe_index].dofs_per_hex), + ExcVectorDoesNotMatch()); + break; + default: + Assert (false, ExcNotImplemented()); + } + // this function really only makes // sense if either a) there are @@ -468,83 +485,50 @@ DoFObjectAccessor<2,DH>::get_dof_indices (std::vector &dof_indices Assert (this->fe_index_is_active (fe_index) || (this->dof_handler->get_fe()[fe_index].dofs_per_cell == - 4*this->dof_handler->get_fe()[fe_index].dofs_per_vertex), + GeometryInfo::vertices_per_cell * + this->dof_handler->get_fe()[fe_index].dofs_per_vertex), ExcInternalError()); - - Assert (static_cast(this->level()) < this->dof_handler->levels.size(), - ExcMessage ("DoFHandler not initialized")); - - const unsigned int dofs_per_vertex = this->dof_handler->get_fe()[fe_index].dofs_per_vertex, - dofs_per_line = this->dof_handler->get_fe()[fe_index].dofs_per_line, - dofs_per_quad = this->dof_handler->get_fe()[fe_index].dofs_per_quad; - std::vector::iterator next = dof_indices.begin(); - for (unsigned int vertex=0; vertex<4; ++vertex) - for (unsigned int d=0; dvertex_dof_index(vertex,d,fe_index); - // now copy dof numbers from the line. for - // lines with the wrong orientation (which - // might occur in 3d), we have already made - // sure that we're ok by picking the correct - // vertices (this happens automatically in - // the vertex() function). however, if the - // line is in wrong orientation, we look at - // it in flipped orientation and we will have - // to adjust the shape function indices that - // we see to correspond to the correct - // (face-local) ordering. - for (unsigned int line=0; line<4; ++line) - for (unsigned int d=0; dline(line)->dof_index(this->dof_handler->get_fe()[fe_index]. - adjust_line_dof_index_for_line_orientation(d, - this->line_orientation(line)),fe_index); - for (unsigned int d=0; ddof_index(d,fe_index); + + // now do the actual work + internal::DoFAccessor::get_dof_indices (*this, dof_indices, fe_index); } -/*------------------------- Functions: DoFObjectAccessor<3,dim> -----------------------*/ -template -inline -DoFObjectAccessor<3,DH>:: -DoFObjectAccessor (const Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data) - : - DoFAccessor<3,DH> (tria, level, index, local_data) -{} - - -template +template inline -TriaIterator > -DoFObjectAccessor<3,DH>::line (const unsigned int i) const +typename internal::DoFHandler::Iterators::line_iterator +DoFAccessor::line (const unsigned int i) const { - TriaIterator > l = TriaObjectAccessor<3,DH::dimension>::line(i); - return TriaIterator > + Assert (structdim > 1, ExcImpossibleInDim(structdim)); + // checking of 'i' happens in + // line_index(i) + + return typename internal::DoFHandler::Iterators::line_iterator ( this->tria, - 0, - l->index(), + 0, // only sub-objects are allowed, which have no level + this->line_index(i), this->dof_handler ); } -template +template inline -TriaIterator > -DoFObjectAccessor<3,DH>::quad (const unsigned int i) const +typename internal::DoFHandler::Iterators::quad_iterator +DoFAccessor::quad (const unsigned int i) const { - Assert (i<6, ExcIndexRange (i, 0, 6)); + Assert (structdim > 2, ExcImpossibleInDim(structdim)); + // checking of 'i' happens in + // quad_index(i) - return TriaIterator > + return typename internal::DoFHandler::Iterators::quad_iterator ( this->tria, - 0, + 0, // only sub-objects are allowed, which have no level this->quad_index (i), this->dof_handler ); @@ -552,119 +536,369 @@ DoFObjectAccessor<3,DH>::quad (const unsigned int i) const -template -inline -void -DoFObjectAccessor<3,DH>::get_dof_indices (std::vector &dof_indices, - const unsigned int fe_index) const -{ - Assert (this->dof_handler != 0, - typename BaseClass::ExcInvalidObject()); - Assert (&this->dof_handler->get_fe() != 0, - typename BaseClass::ExcInvalidObject()); - Assert (dof_indices.size() == (8*this->dof_handler->get_fe()[fe_index].dofs_per_vertex + - 12*this->dof_handler->get_fe()[fe_index].dofs_per_line + - 6*this->dof_handler->get_fe()[fe_index].dofs_per_quad + - this->dof_handler->get_fe()[fe_index].dofs_per_hex), - typename BaseClass::ExcVectorDoesNotMatch()); - Assert (static_cast(this->level()) < this->dof_handler->levels.size(), - ExcMessage ("DoFHandler not initialized")); - // this function really only makes - // sense if either a) there are - // degrees of freedom defined on - // the present object, or b) the - // object is non-active objects but - // all degrees of freedom are - // located on vertices, since - // otherwise there are degrees of - // freedom on sub-objects which are - // not allocated for this - // non-active thing - Assert (this->fe_index_is_active (fe_index) - || - (this->dof_handler->get_fe()[fe_index].dofs_per_cell == - 8*this->dof_handler->get_fe()[fe_index].dofs_per_vertex), - ExcInternalError()); - Assert (static_cast(this->level()) < this->dof_handler->levels.size(), - ExcMessage ("DoFHandler not initialized")); - - const unsigned int dofs_per_vertex = this->dof_handler->get_fe()[fe_index].dofs_per_vertex, - dofs_per_line = this->dof_handler->get_fe()[fe_index].dofs_per_line, - dofs_per_quad = this->dof_handler->get_fe()[fe_index].dofs_per_quad, - dofs_per_hex = this->dof_handler->get_fe()[fe_index].dofs_per_hex; - std::vector::iterator next = dof_indices.begin(); - for (unsigned int vertex=0; vertex<8; ++vertex) - for (unsigned int d=0; dvertex_dof_index(vertex,d,fe_index); - // now copy dof numbers from the line. for - // lines with the wrong orientation, we have - // already made sure that we're ok by picking - // the correct vertices (this happens - // automatically in the vertex() - // function). however, if the line is in - // wrong orientation, we look at it in - // flipped orientation and we will have to - // adjust the shape function indices that we - // see to correspond to the correct - // (cell-local) ordering. - for (unsigned int line=0; line<12; ++line) - for (unsigned int d=0; dline(line)->dof_index(this->dof_handler->get_fe()[fe_index]. - adjust_line_dof_index_for_line_orientation(d, - this->line_orientation(line)),fe_index); - // now copy dof numbers from the face. for - // faces with the wrong orientation, we - // have already made sure that we're ok by - // picking the correct lines and vertices - // (this happens automatically in the - // line() and vertex() functions). however, - // if the face is in wrong orientation, we - // look at it in flipped orientation and we - // will have to adjust the shape function - // indices that we see to correspond to the - // correct (cell-local) ordering. The same - // applies, if the face_rotation or - // face_orientation is non-standard - for (unsigned int quad=0; quad<6; ++quad) - for (unsigned int d=0; dquad(quad)->dof_index(this->dof_handler->get_fe()[fe_index]. - adjust_quad_dof_index_for_face_orientation(d, - this->face_orientation(quad), - this->face_flip(quad), - this->face_rotation(quad)),fe_index); - for (unsigned int d=0; ddof_index(d,fe_index); -} +/*------------------------- Functions: DoFCellAccessor -----------------------*/ +/** + * A class with the same purpose as the similarly named class of the + * Triangulation class. See there for more information. + */ +template +struct DoFCellAccessor::Implementation +{ + /** + * Implement the updating of the + * cache. Currently not + * implemented for hp::DoFHandler + * objects. + */ + template + static + void + update_cell_dof_indices_cache (const DoFCellAccessor > &accessor) + { + // check as in documentation that + // cell is either active, or dofs + // are only in vertices. otherwise + // simply don't update the cache at + // all. the get_dof_indices + // function will then make sure we + // don't access the invalid data + if (accessor.has_children() + && + (accessor.get_fe().dofs_per_cell != + accessor.get_fe().dofs_per_vertex * GeometryInfo<1>::vertices_per_cell)) + return; + + const unsigned int dofs_per_vertex = accessor.get_fe().dofs_per_vertex, + dofs_per_line = accessor.get_fe().dofs_per_line, + dofs_per_cell = accessor.get_fe().dofs_per_cell; + + // make sure the cache is at least + // as big as we need it when + // writing to the last element of + // this cell + Assert (accessor.present_index * dofs_per_cell + dofs_per_cell + <= + accessor.dof_handler->levels[accessor.present_level] + ->cell_dof_indices_cache.size(), + ExcInternalError()); + + std::vector::iterator next + = (accessor.dof_handler->levels[accessor.present_level] + ->cell_dof_indices_cache.begin() + accessor.present_index * dofs_per_cell); + + for (unsigned int vertex=0; vertex<2; ++vertex) + for (unsigned int d=0; d + static + void + update_cell_dof_indices_cache (const DoFCellAccessor > &accessor) + { + // check as in documentation that + // cell is either active, or dofs + // are only in vertices. otherwise + // simply don't update the cache at + // all. the get_dof_indices + // function will then make sure we + // don't access the invalid data + if (accessor.has_children() + && + (accessor.get_fe().dofs_per_cell != + accessor.get_fe().dofs_per_vertex * GeometryInfo<2>::vertices_per_cell)) + return; + + const unsigned int dofs_per_vertex = accessor.get_fe().dofs_per_vertex, + dofs_per_line = accessor.get_fe().dofs_per_line, + dofs_per_quad = accessor.get_fe().dofs_per_quad, + dofs_per_cell = accessor.get_fe().dofs_per_cell; + + // make sure the cache is at least + // as big as we need it when + // writing to the last element of + // this cell + Assert (accessor.present_index * dofs_per_cell + dofs_per_cell + <= + accessor.dof_handler->levels[accessor.present_level] + ->cell_dof_indices_cache.size(), + ExcInternalError()); + + std::vector::iterator next + = (accessor.dof_handler->levels[accessor.present_level] + ->cell_dof_indices_cache.begin() + accessor.present_index * dofs_per_cell); + + for (unsigned int vertex=0; vertex<4; ++vertex) + for (unsigned int d=0; ddof_index(d); + for (unsigned int d=0; d + static + void + update_cell_dof_indices_cache (const DoFCellAccessor > &accessor) + { + // check as in documentation that + // cell is either active, or dofs + // are only in vertices. otherwise + // simply don't update the cache at + // all. the get_dof_indices + // function will then make sure we + // don't access the invalid data + if (accessor.has_children() + && + (accessor.get_fe().dofs_per_cell != + accessor.get_fe().dofs_per_vertex * GeometryInfo<3>::vertices_per_cell)) + return; + + const unsigned int dofs_per_vertex = accessor.get_fe().dofs_per_vertex, + dofs_per_line = accessor.get_fe().dofs_per_line, + dofs_per_quad = accessor.get_fe().dofs_per_quad, + dofs_per_hex = accessor.get_fe().dofs_per_hex, + dofs_per_cell = accessor.get_fe().dofs_per_cell; + + // make sure the cache is at least + // as big as we need it when + // writing to the last element of + // this cell + Assert (accessor.present_index * dofs_per_cell + dofs_per_cell + <= + accessor.dof_handler->levels[accessor.present_level] + ->cell_dof_indices_cache.size(), + ExcInternalError()); + + std::vector::iterator next + = (accessor.dof_handler->levels[accessor.present_level] + ->cell_dof_indices_cache.begin() + accessor.present_index * dofs_per_cell); + + for (unsigned int vertex=0; vertex<8; ++vertex) + for (unsigned int d=0; ddof_index(accessor.dof_handler->get_fe(). + adjust_line_dof_index_for_line_orientation(d, + accessor.line_orientation(line))); + // now copy dof numbers from the face. for + // faces with the wrong orientation, we + // have already made sure that we're ok by + // picking the correct lines and vertices + // (this happens automatically in the + // line() and vertex() functions). however, + // if the face is in wrong orientation, we + // look at it in flipped orientation and we + // will have to adjust the shape function + // indices that we see to correspond to the + // correct (cell-local) ordering. The same + // applies, if the face_rotation or + // face_orientation is non-standard + for (unsigned int quad=0; quad<6; ++quad) + for (unsigned int d=0; ddof_index(accessor.dof_handler->get_fe(). + adjust_quad_dof_index_for_face_orientation(d, + accessor.face_orientation(quad), + accessor.face_flip(quad), + accessor.face_rotation(quad))); + for (unsigned int d=0; d + static + void + update_cell_dof_indices_cache (const DoFCellAccessor > &) + { +//TODO[WB]: should implement a dof indices cache for hp as well + + // not implemented, but should also + // not be called + Assert (false, ExcNotImplemented()); + } + + + /** + * A function that collects the + * global indices of degrees of + * freedom. This function works + * for ::DoFHandler and all + * template arguments and copies + * the data out of the cache that + * we hold for each cell. + */ + template + static + void + get_dof_indices (const DoFCellAccessor > &accessor, + std::vector &dof_indices) + { + Assert (dof_indices.size() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + + // check as in documentation that + // cell is either active, or dofs + // are only in vertices + Assert (!accessor.has_children() + || + (accessor.get_fe().dofs_per_cell == + accessor.get_fe().dofs_per_vertex * GeometryInfo::vertices_per_cell), + ExcMessage ("Cell must either be active, or all DoFs must be in vertices")); + + unsigned int *cache = &accessor.dof_handler->levels[accessor.level()] + ->cell_dof_indices_cache[accessor.present_index * accessor.get_fe().dofs_per_cell]; + for (unsigned int i=0; i + static + void + get_dof_indices (const DoFCellAccessor > &accessor, + std::vector &dof_indices) + { + // no caching for hp::DoFHandler implemented + accessor.DoFAccessor >::get_dof_indices (dof_indices, + accessor.active_fe_index()); + } + + + /** + * Do what the active_fe_index + * function in the parent class + * is supposed to do. + */ + template + static + unsigned int + active_fe_index (const DoFCellAccessor > &) + { + // ::DoFHandler only supports a + // single active fe with index + // zero + return 0; + } + + + + template + static + unsigned int + active_fe_index (const DoFCellAccessor > &accessor) + { + Assert (static_cast(accessor.level()) < accessor.dof_handler->levels.size(), + ExcMessage ("DoFHandler not initialized")); + Assert (static_cast::size_type>(accessor.present_index) < + accessor.dof_handler->levels[accessor.level()]->active_fe_indices.size (), + ExcIndexRange (accessor.present_index, 0, + accessor.dof_handler->levels[accessor.level()]->active_fe_indices.size ())); + return accessor.dof_handler->levels[accessor.level()] + ->active_fe_indices[accessor.present_index]; + } + + + + /** + * Do what the + * set_active_fe_index function + * in the parent class is + * supposed to do. + */ + template + static + void + set_active_fe_index (const DoFCellAccessor > &, + const unsigned int i) + { + // ::DoFHandler only supports a + // single active fe with index + // zero + Assert (i == 0, typename BaseClass::ExcInvalidObject()); + } + + + + template + static + void + set_active_fe_index (DoFCellAccessor > &accessor, + const unsigned int i) + { + Assert (accessor.dof_handler != 0, + typename BaseClass::ExcInvalidObject()); + Assert (static_cast(accessor.level()) < + accessor.dof_handler->levels.size(), + ExcMessage ("DoFHandler not initialized")); + Assert (static_cast::size_type>(accessor.present_index) < + accessor.dof_handler->levels[accessor.level()]->active_fe_indices.size (), + ExcIndexRange (accessor.present_index, 0, + accessor.dof_handler->levels[accessor.level()]->active_fe_indices.size ())); + accessor.dof_handler->levels[accessor.level()] + ->active_fe_indices[accessor.present_index] = i; + } +}; -/*------------------------- Functions: DoFCellAccessor -----------------------*/ template inline DoFCellAccessor:: -DoFCellAccessor (const Triangulation *tria, +DoFCellAccessor (const Triangulation *tria, const int level, const int index, const AccessorData *local_data) : - DoFObjectAccessor (tria,level,index,local_data) + DoFAccessor (tria,level,index,local_data) {} template inline -TriaIterator > +typename internal::DoFHandler::Iterators::cell_iterator DoFCellAccessor::neighbor (const unsigned int i) const { - TriaIterator > q (this->tria, - this->neighbor_level (i), - this->neighbor_index (i), - this->dof_handler); + typename internal::DoFHandler::Iterators::cell_iterator + q (this->tria, + this->neighbor_level (i), + this->neighbor_index (i), + this->dof_handler); #ifdef DEBUG if (q.state() != IteratorState::past_the_end) @@ -676,13 +910,14 @@ DoFCellAccessor::neighbor (const unsigned int i) const template inline -TriaIterator > +typename internal::DoFHandler::Iterators::cell_iterator DoFCellAccessor::child (const unsigned int i) const { - TriaIterator > q (this->tria, - this->level()+1, - this->child_index (i), - this->dof_handler); + typename internal::DoFHandler::Iterators::cell_iterator + q (this->tria, + this->level()+1, + this->child_index (i), + this->dof_handler); #ifdef DEBUG if (q.state() != IteratorState::past_the_end) @@ -693,159 +928,58 @@ DoFCellAccessor::child (const unsigned int i) const -template <> -inline -TriaIterator<1, DoFObjectAccessor<0,hp::DoFHandler<1> > > -DoFCellAccessor >::face (const unsigned int) const -{ - Assert (false, ExcImpossibleInDim(1)); - return TriaIterator<1, DoFObjectAccessor<0,hp::DoFHandler<1> > >(); -} - - - -template <> -inline -TriaIterator<2, DoFObjectAccessor<1,hp::DoFHandler<2> > > -DoFCellAccessor >::face (const unsigned int i) const -{ - return this->line(i); -} - - - -template <> -inline -TriaIterator<3, DoFObjectAccessor<2, hp::DoFHandler<3> > > -DoFCellAccessor >::face (const unsigned int i) const -{ - return this->quad(i); -} - - - -template <> -inline -void -DoFCellAccessor >:: -get_dof_indices (std::vector &dof_indices) const -{ - Assert (dof_indices.size() == this->get_fe().dofs_per_cell, - BaseClass::ExcVectorDoesNotMatch()); - - // check as in documentation that - // cell is either active, or dofs - // are only in vertices - Assert (!this->has_children() - || - (this->get_fe().dofs_per_cell == - this->get_fe().dofs_per_vertex * GeometryInfo<1>::vertices_per_cell), - ExcMessage ("Cell must either be active, or all DoFs must be in vertices")); - - unsigned int *cache = &this->dof_handler->levels[this->level()] - ->cell_dof_indices_cache[this->present_index * this->get_fe().dofs_per_cell]; - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i, ++cache) - dof_indices[i] = *cache; -} - - - -template <> -inline -void -DoFCellAccessor >:: -get_dof_indices (std::vector &dof_indices) const -{ - // no caching for hp::DoFHandler implemented - DoFObjectAccessor >::get_dof_indices (dof_indices, - this->active_fe_index()); -} - - -template <> -inline -void -DoFCellAccessor >:: -get_dof_indices (std::vector &dof_indices) const -{ - Assert (dof_indices.size() == this->get_fe().dofs_per_cell, - BaseClass::ExcVectorDoesNotMatch()); - - // check as in documentation that - // cell is either active, or dofs - // are only in vertices - Assert (!this->has_children() - || - (this->get_fe().dofs_per_cell == - this->get_fe().dofs_per_vertex * GeometryInfo<2>::vertices_per_cell), - ExcMessage ("Cell must either be active, or all DoFs must be in vertices")); - - unsigned int *cache = &this->dof_handler->levels[this->level()] - ->cell_dof_indices_cache[this->present_index * this->get_fe().dofs_per_cell]; - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i, ++cache) - dof_indices[i] = *cache; -} - - - -template <> -inline -void -DoFCellAccessor >:: -get_dof_indices (std::vector &dof_indices) const -{ - // no caching for hp::DoFHandler implemented - DoFObjectAccessor >::get_dof_indices (dof_indices, - this->active_fe_index()); -} - - - -template <> -inline -void -DoFCellAccessor >:: -get_dof_indices (std::vector &dof_indices) const +template +typename internal::DoFHandler::Iterators::face_iterator +DoFCellAccessor::face (const unsigned int i) const { - Assert (dof_indices.size() == this->get_fe().dofs_per_cell, - BaseClass::ExcVectorDoesNotMatch()); + Assert (i::faces_per_cell, ExcIndexRange (i, 0, GeometryInfo::faces_per_cell)); + Assert (static_cast(this->level()) < this->dof_handler->levels.size(), + ExcMessage ("DoFHandler not initialized")); - // check as in documentation that - // cell is either active, or dofs - // are only in vertices - Assert (!this->has_children() - || - (this->get_fe().dofs_per_cell == - this->get_fe().dofs_per_vertex * GeometryInfo<3>::vertices_per_cell), - ExcMessage ("Cell must either be active, or all DoFs must be in vertices")); - - unsigned int *cache = &this->dof_handler->levels[this->level()] - ->cell_dof_indices_cache[this->present_index * this->get_fe().dofs_per_cell]; - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i, ++cache) - dof_indices[i] = *cache; + const unsigned int dim = DH::dimension; + Assert (dim > 1, ExcImpossibleInDim(1)); + + switch (dim) + { + case 2: + return typename internal::DoFHandler::Iterators::face_iterator + (this->tria, + 0, + this->line_index (i), + this->dof_handler); + + case 3: + return typename internal::DoFHandler::Iterators::face_iterator + (this->tria, + 0, + this->quad_index (i), + this->dof_handler); + + default: + Assert (false, ExcNotImplemented()); + return typename internal::DoFHandler::Iterators::face_iterator(); + } } -template <> +template inline void -DoFCellAccessor >:: +DoFCellAccessor:: get_dof_indices (std::vector &dof_indices) const { - // no caching for hp::DoFHandler implemented - DoFObjectAccessor >::get_dof_indices (dof_indices, - this->active_fe_index()); + Implementation::get_dof_indices (*this, dof_indices); } template inline -const FiniteElement & +const FiniteElement & DoFCellAccessor::get_fe () const { - return this->dof_handler->get_fe(); + return internal::DoFAccessor::get_fe (this->dof_handler->get_fe(), active_fe_index()); } @@ -855,7 +989,7 @@ inline unsigned int DoFCellAccessor::active_fe_index () const { - return 0; + return Implementation::active_fe_index (*this); } @@ -865,153 +999,11 @@ inline void DoFCellAccessor::set_active_fe_index (const unsigned int i) { - Assert (i == 0, typename BaseClass::BaseClass::ExcInvalidObject()); -} - - - - - -template <> -inline -const FiniteElement<1> & -DoFCellAccessor >::get_fe () const -{ - return this->dof_handler->get_fe()[active_fe_index ()]; -} - - - -template <> -inline -const FiniteElement<2> & -DoFCellAccessor >::get_fe () const -{ - return this->dof_handler->get_fe()[active_fe_index ()]; + Implementation::set_active_fe_index (*this, i); } -template <> -inline -const FiniteElement<3> & -DoFCellAccessor >::get_fe () const -{ - return this->dof_handler->get_fe()[active_fe_index ()]; -} - - - -template <> -inline -unsigned int -DoFCellAccessor >::active_fe_index () const -{ - Assert (static_cast(this->level()) < this->dof_handler->levels.size(), - ExcMessage ("DoFHandler not initialized")); - Assert (static_cast::size_type>(this->present_index) < - this->dof_handler->levels[this->level()]->active_fe_indices.size (), - ExcIndexRange (this->present_index, 0, - this->dof_handler->levels[this->level()]->active_fe_indices.size ())); - return this->dof_handler->levels[this->level()] - ->active_fe_indices[this->present_index]; -} - - - -template <> -inline -unsigned int -DoFCellAccessor >::active_fe_index () const -{ - Assert (static_cast(this->level()) < this->dof_handler->levels.size(), - ExcMessage ("DoFHandler not initialized")); - Assert (static_cast::size_type>(this->present_index) < - this->dof_handler->levels[this->level()]->active_fe_indices.size (), - ExcIndexRange (this->present_index, 0, - this->dof_handler->levels[this->level()]->active_fe_indices.size ())); - return this->dof_handler->levels[this->level()] - ->active_fe_indices[this->present_index]; -} - - - -template <> -inline -unsigned int -DoFCellAccessor >::active_fe_index () const -{ - Assert (static_cast(this->level()) < this->dof_handler->levels.size(), - ExcMessage ("DoFHandler not initialized")); - Assert (static_cast::size_type>(this->present_index) < - this->dof_handler->levels[this->level()]->active_fe_indices.size (), - ExcIndexRange (this->present_index, 0, - this->dof_handler->levels[this->level()]->active_fe_indices.size ())); - return this->dof_handler->levels[this->level()] - ->active_fe_indices[this->present_index]; -} - - - -template <> -inline -void -DoFCellAccessor >::set_active_fe_index (const unsigned int i) -{ - Assert (this->dof_handler != 0, - BaseClass::ExcInvalidObject()); - Assert (static_cast(this->level()) < - this->dof_handler->levels.size(), - ExcMessage ("DoFHandler not initialized")); - Assert (static_cast::size_type>(this->present_index) < - this->dof_handler->levels[this->level()]->active_fe_indices.size (), - ExcIndexRange (this->present_index, 0, - this->dof_handler->levels[this->level()]->active_fe_indices.size ())); - this->dof_handler->levels[this->level()] - ->active_fe_indices[this->present_index] = i; -} - - - -template <> -inline -void -DoFCellAccessor >::set_active_fe_index (const unsigned int i) -{ - Assert (this->dof_handler != 0, - BaseClass::ExcInvalidObject()); - Assert (static_cast(this->level()) < - this->dof_handler->levels.size(), - ExcMessage ("DoFHandler not initialized")); - Assert (static_cast::size_type>(this->present_index) < - this->dof_handler->levels[this->level()]->active_fe_indices.size (), - ExcIndexRange (this->present_index, 0, - this->dof_handler->levels[this->level()]->active_fe_indices.size ())); - this->dof_handler->levels[this->level()] - ->active_fe_indices[this->present_index] = i; -} - - - -template <> -inline -void -DoFCellAccessor >::set_active_fe_index (const unsigned int i) -{ - Assert (this->dof_handler != 0, - BaseClass::ExcInvalidObject()); - Assert (static_cast(this->level()) < - this->dof_handler->levels.size(), - ExcMessage ("DoFHandler not initialized")); - Assert (static_cast::size_type>(this->present_index) < - this->dof_handler->levels[this->level()]->active_fe_indices.size (), - ExcIndexRange (this->present_index, 0, - this->dof_handler->levels[this->level()]->active_fe_indices.size ())); - this->dof_handler->levels[this->level()] - ->active_fe_indices[this->present_index] = i; -} - - template template diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index b6a757bf3a..176092f490 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -48,7 +48,7 @@ class BlockIndices; * $x_{i_1} = \sum_{j=2}^M a_{i_j} x_{i_j}$. In the context of adaptive finite * elements, such constraints appear most frequently as "hanging nodes". For * example, when using Q1 and Q2 elements (i.e. using - * FE_Q<dim>(1) and FE_Q<dim>(2)) on the two + * FE_Q<dim,spacedim>(1) and FE_Q<dim,spacedim>(2)) on the two * marked cells of the mesh * * @image html hp-refinement-simple.png diff --git a/deal.II/deal.II/include/dofs/dof_faces.h b/deal.II/deal.II/include/dofs/dof_faces.h index 658ef72767..33e51963b5 100644 --- a/deal.II/deal.II/include/dofs/dof_faces.h +++ b/deal.II/deal.II/include/dofs/dof_faces.h @@ -21,8 +21,8 @@ DEAL_II_NAMESPACE_OPEN -template class DoFHandler; -template class MGDoFHandler; +template class DoFHandler; +template class MGDoFHandler; namespace internal { diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index 1f9ec0668f..cab7749455 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -62,6 +62,13 @@ namespace internal * presented by the ++ and -- operators is the same as that * for the corresponding triangulation iterators. * + * The spacedim parameter has to be used if one wants to + * solve problems in the boundary element method formulation or in an + * equivalent one, as it is explained in the Triangulation class. If + * not specified, this parameter takes the default value =dim + * so that this class can be used to solve problems in the finite + * element method formulation. + * * *

Distribution of indices for degrees of freedom

* @@ -104,10 +111,10 @@ namespace internal * @ingroup dofs * @author Wolfgang Bangerth, 1998 */ -template +template class DoFHandler : public Subscriptor { - typedef internal::DoFHandler::Iterators > IteratorSelector; + typedef internal::DoFHandler::Iterators > IteratorSelector; public: typedef typename IteratorSelector::raw_line_iterator raw_line_iterator; typedef typename IteratorSelector::line_iterator line_iterator; @@ -133,14 +140,20 @@ class DoFHandler : public Subscriptor * Alias the @p FunctionMap type * declared elsewhere. */ - typedef typename dealii::FunctionMap::type FunctionMap; + typedef typename dealii::FunctionMap::type FunctionMap; /** * Make the dimension available * in function templates. */ static const unsigned int dimension = dim; - + + /** + * Make the space dimension available + * in function templates. + */ + static const unsigned int space_dimension = spacedim; + /** * When the arrays holding the * DoF indices are set up, but @@ -184,8 +197,8 @@ class DoFHandler : public Subscriptor * Constructor. Take @p tria as the * triangulation to work on. */ - DoFHandler (const Triangulation &tria); - + DoFHandler ( const Triangulation &tria); + /** * Destructor. */ @@ -222,7 +235,7 @@ class DoFHandler : public Subscriptor * releases the lock of this * object to the finite element. */ - virtual void distribute_dofs (const FiniteElement &fe, + virtual void distribute_dofs (const FiniteElement &fe, const unsigned int offset = 0); /** @@ -898,14 +911,14 @@ class DoFHandler : public Subscriptor * the selected finite element * object. */ - const FiniteElement & get_fe () const; + const FiniteElement & get_fe () const; /** * Return a constant reference to * the triangulation underlying * this object. */ - const Triangulation & get_tria () const; + const Triangulation & get_tria () const; /** * Determine an estimate for the @@ -954,6 +967,25 @@ class DoFHandler : public Subscriptor int, << "The given list of new dof indices is not consecutive: " << "the index " << arg1 << " does not exist."); + /** + * Exception + */ + DeclException1 (ExcInvalidLevel, + int, + << "The given level " << arg1 + << " is not in the valid range!"); + /** + * Exception + */ + DeclException0 (ExcFacesHaveNoLevel); + /** + * The triangulation level you + * accessed is empty. + */ + DeclException1 (ExcEmptyLevel, + int, + << "You tried to do something on level " << arg1 + << ", but this level is empty."); protected: @@ -961,7 +993,7 @@ class DoFHandler : public Subscriptor * Address of the triangulation to * work on. */ - SmartPointer > tria; + SmartPointer > tria; /** * Store a pointer to the finite element @@ -975,7 +1007,7 @@ class DoFHandler : public Subscriptor * function (this clears all data of * this object as well, though). */ - SmartPointer > selected_fe; + SmartPointer > selected_fe; private: @@ -1001,38 +1033,10 @@ class DoFHandler : public Subscriptor */ DoFHandler & operator = (const DoFHandler &); - /** - * Reserve enough space in the - * levels[] objects to store the - * numbers of the degrees of freedom - * needed for the given element. The - * given element is that one which - * was selected when calling - * @p distribute_dofs the last time. - */ - void reserve_space (); - /** * Free all used memory. */ void clear_space (); - - /** - * Distribute dofs on the given cell, - * with new dofs starting with index - * @p next_free_dof. Return the next - * unused index number. The finite - * element used is the one given to - * @p distribute_dofs, which is copied - * to @p selected_fe. - * - * This function is excluded from the - * @p distribute_dofs function since - * it can not be implemented dimension - * independent. - */ - unsigned int distribute_dofs_on_cell (active_cell_iterator &cell, - unsigned int next_free_dof); /** * Set the @p local_index-th @@ -1199,21 +1203,28 @@ class DoFHandler : public Subscriptor */ std::vector vertex_dofs; + /** + * Forward declaration of a class + * into which we put significant + * parts of the implementation. + * + * See the .cc file for more + * information. + */ + struct Implementation; /* * Make accessor objects friends. */ template friend class DoFAccessor; - /* - * Make accessor objects friends. - */ - template friend class DoFObjectAccessor; - /* * Same with cell accessor */ template friend class DoFCellAccessor; + + + friend class DoFHandler::Implementation; }; @@ -1226,153 +1237,50 @@ class DoFHandler : public Subscriptor template <> unsigned int DoFHandler<1>::n_boundary_dofs () const; template <> unsigned int DoFHandler<1>::n_boundary_dofs (const FunctionMap &) const; template <> unsigned int DoFHandler<1>::n_boundary_dofs (const std::set &) const; -template <> unsigned int DoFHandler<1>::max_couplings_between_dofs () const; -template <> unsigned int DoFHandler<1>::max_couplings_between_boundary_dofs () const; -template <> unsigned int DoFHandler<2>::max_couplings_between_dofs () const; -template <> unsigned int DoFHandler<2>::max_couplings_between_boundary_dofs () const; -template <> unsigned int DoFHandler<3>::max_couplings_between_dofs () const; -template <> unsigned int DoFHandler<3>::max_couplings_between_boundary_dofs () const; - -template <> DoFHandler<1>::raw_cell_iterator DoFHandler<1>::begin_raw (const unsigned int level) const; -template <> DoFHandler<1>::cell_iterator DoFHandler<1>::begin (const unsigned int level) const; -template <> DoFHandler<1>::active_cell_iterator DoFHandler<1>::begin_active (const unsigned int level) const; -template <> DoFHandler<1>::raw_cell_iterator DoFHandler<1>::end () const; -template <> DoFHandler<1>::raw_cell_iterator DoFHandler<1>::last_raw () const; -template <> DoFHandler<1>::raw_cell_iterator DoFHandler<1>::last_raw (const unsigned int level) const; -template <> DoFHandler<1>::cell_iterator DoFHandler<1>::last () const; -template <> DoFHandler<1>::cell_iterator DoFHandler<1>::last (const unsigned int level) const; -template <> DoFHandler<1>::active_cell_iterator DoFHandler<1>::last_active () const; -template <> DoFHandler<1>::active_cell_iterator DoFHandler<1>::last_active (const unsigned int level) const; -template <> DoFHandler<1>::raw_face_iterator DoFHandler<1>::begin_raw_face () const; -template <> DoFHandler<1>::face_iterator DoFHandler<1>::begin_face () const; -template <> DoFHandler<1>::active_face_iterator DoFHandler<1>::begin_active_face () const; -template <> DoFHandler<1>::raw_face_iterator DoFHandler<1>::end_face () const; -template <> DoFHandler<1>::raw_face_iterator DoFHandler<1>::last_raw_face () const; -template <> DoFHandler<1>::face_iterator DoFHandler<1>::last_face () const; -template <> DoFHandler<1>::active_face_iterator DoFHandler<1>::last_active_face () const; -template <> DoFHandler<1>::raw_quad_iterator DoFHandler<1>::begin_raw_quad (const unsigned int) const; -template <> DoFHandler<1>::quad_iterator DoFHandler<1>::begin_quad (const unsigned int) const; -template <> DoFHandler<1>::active_quad_iterator DoFHandler<1>::begin_active_quad (const unsigned int) const; -template <> DoFHandler<1>::raw_quad_iterator DoFHandler<1>::end_quad () const; -template <> DoFHandler<1>::raw_quad_iterator DoFHandler<1>::last_raw_quad (const unsigned int) const; -template <> DoFHandler<1>::quad_iterator DoFHandler<1>::last_quad (const unsigned int) const; -template <> DoFHandler<1>::active_quad_iterator DoFHandler<1>::last_active_quad (const unsigned int) const; -template <> DoFHandler<1>::raw_quad_iterator DoFHandler<1>::last_raw_quad () const; -template <> DoFHandler<1>::quad_iterator DoFHandler<1>::last_quad () const; -template <> DoFHandler<1>::active_quad_iterator DoFHandler<1>::last_active_quad () const; -template <> DoFHandler<1>::raw_hex_iterator DoFHandler<1>::begin_raw_hex (const unsigned int) const; -template <> DoFHandler<1>::hex_iterator DoFHandler<1>::begin_hex (const unsigned int) const; -template <> DoFHandler<1>::active_hex_iterator DoFHandler<1>::begin_active_hex (const unsigned int) const; -template <> DoFHandler<1>::raw_hex_iterator DoFHandler<1>::end_hex () const; -template <> DoFHandler<1>::raw_hex_iterator DoFHandler<1>::last_raw_hex (const unsigned int) const; -template <> DoFHandler<1>::raw_hex_iterator DoFHandler<1>::last_raw_hex () const; -template <> DoFHandler<1>::hex_iterator DoFHandler<1>::last_hex (const unsigned int) const; -template <> DoFHandler<1>::hex_iterator DoFHandler<1>::last_hex () const; -template <> DoFHandler<1>::active_hex_iterator DoFHandler<1>::last_active_hex (const unsigned int) const; -template <> DoFHandler<1>::active_hex_iterator DoFHandler<1>::last_active_hex () const; -template <> DoFHandler<2>::raw_cell_iterator DoFHandler<2>::begin_raw (const unsigned int level) const; -template <> DoFHandler<2>::cell_iterator DoFHandler<2>::begin (const unsigned int level) const; -template <> DoFHandler<2>::active_cell_iterator DoFHandler<2>::begin_active (const unsigned int level) const; -template <> DoFHandler<2>::raw_cell_iterator DoFHandler<2>::end () const; -template <> DoFHandler<2>::raw_cell_iterator DoFHandler<2>::last_raw () const; -template <> DoFHandler<2>::raw_cell_iterator DoFHandler<2>::last_raw (const unsigned int level) const; -template <> DoFHandler<2>::cell_iterator DoFHandler<2>::last () const; -template <> DoFHandler<2>::cell_iterator DoFHandler<2>::last (const unsigned int level) const; -template <> DoFHandler<2>::active_cell_iterator DoFHandler<2>::last_active () const; -template <> DoFHandler<2>::active_cell_iterator DoFHandler<2>::last_active (const unsigned int level) const; -template <> DoFHandler<2>::raw_face_iterator DoFHandler<2>::begin_raw_face () const; -template <> DoFHandler<2>::face_iterator DoFHandler<2>::begin_face () const; -template <> DoFHandler<2>::active_face_iterator DoFHandler<2>::begin_active_face () const; -template <> DoFHandler<2>::raw_face_iterator DoFHandler<2>::end_face () const; -template <> DoFHandler<2>::raw_face_iterator DoFHandler<2>::last_raw_face () const; -template <> DoFHandler<2>::face_iterator DoFHandler<2>::last_face () const; -template <> DoFHandler<2>::active_face_iterator DoFHandler<2>::last_active_face () const; -template <> DoFHandler<2>::raw_hex_iterator DoFHandler<2>::begin_raw_hex (const unsigned int) const; -template <> DoFHandler<2>::hex_iterator DoFHandler<2>::begin_hex (const unsigned int) const; -template <> DoFHandler<2>::active_hex_iterator DoFHandler<2>::begin_active_hex (const unsigned int) const; -template <> DoFHandler<2>::raw_hex_iterator DoFHandler<2>::end_hex () const; -template <> DoFHandler<2>::raw_hex_iterator DoFHandler<2>::last_raw_hex (const unsigned int) const; -template <> DoFHandler<2>::raw_hex_iterator DoFHandler<2>::last_raw_hex () const; -template <> DoFHandler<2>::hex_iterator DoFHandler<2>::last_hex (const unsigned int) const; -template <> DoFHandler<2>::hex_iterator DoFHandler<2>::last_hex () const; -template <> DoFHandler<2>::active_hex_iterator DoFHandler<2>::last_active_hex (const unsigned int) const; -template <> DoFHandler<2>::active_hex_iterator DoFHandler<2>::last_active_hex () const; -template <> DoFHandler<3>::raw_cell_iterator DoFHandler<3>::begin_raw (const unsigned int level) const; -template <> DoFHandler<3>::cell_iterator DoFHandler<3>::begin (const unsigned int level) const; -template <> DoFHandler<3>::active_cell_iterator DoFHandler<3>::begin_active (const unsigned int level) const; -template <> DoFHandler<3>::raw_cell_iterator DoFHandler<3>::end () const; -template <> DoFHandler<3>::raw_cell_iterator DoFHandler<3>::last_raw () const; -template <> DoFHandler<3>::raw_cell_iterator DoFHandler<3>::last_raw (const unsigned int level) const; -template <> DoFHandler<3>::cell_iterator DoFHandler<3>::last () const; -template <> DoFHandler<3>::cell_iterator DoFHandler<3>::last (const unsigned int level) const; -template <> DoFHandler<3>::active_cell_iterator DoFHandler<3>::last_active () const; -template <> DoFHandler<3>::active_cell_iterator DoFHandler<3>::last_active (const unsigned int level) const; -template <> DoFHandler<3>::raw_face_iterator DoFHandler<3>::begin_raw_face () const; -template <> DoFHandler<3>::face_iterator DoFHandler<3>::begin_face () const; -template <> DoFHandler<3>::active_face_iterator DoFHandler<3>::begin_active_face () const; -template <> DoFHandler<3>::raw_face_iterator DoFHandler<3>::end_face () const; -template <> DoFHandler<3>::raw_face_iterator DoFHandler<3>::last_raw_face () const; -template <> DoFHandler<3>::face_iterator DoFHandler<3>::last_face () const; -template <> DoFHandler<3>::active_face_iterator DoFHandler<3>::last_active_face () const; - -template <> -unsigned int DoFHandler<1>::distribute_dofs_on_cell (active_cell_iterator &cell, - unsigned int next_free_dof); -template <> -unsigned int DoFHandler<2>::distribute_dofs_on_cell (active_cell_iterator &cell, - unsigned int next_free_dof); -template <> -unsigned int DoFHandler<3>::distribute_dofs_on_cell (active_cell_iterator &cell, - unsigned int next_free_dof); -template <> void DoFHandler<1>::renumber_dofs (const std::vector &new_numbers); -template <> void DoFHandler<2>::renumber_dofs (const std::vector &new_numbers); -template <> void DoFHandler<3>::renumber_dofs (const std::vector &new_numbers); -template <> void DoFHandler<1>::reserve_space (); -template <> void DoFHandler<2>::reserve_space (); -template <> void DoFHandler<3>::reserve_space (); + /* ----------------------- Inline functions ---------------------------------- */ -template +template inline unsigned int -DoFHandler::n_dofs () const +DoFHandler::n_dofs () const { return used_dofs; } -template +template inline -const FiniteElement & -DoFHandler::get_fe () const +const FiniteElement & +DoFHandler::get_fe () const { Assert(selected_fe!=0, ExcNoFESelected()); return *selected_fe; } -template +template inline -const Triangulation & -DoFHandler::get_tria () const +const Triangulation & +DoFHandler::get_tria () const { return *tria; } -template +template inline unsigned int -DoFHandler:: +DoFHandler:: get_vertex_dof_index (const unsigned int vertex_index, const unsigned int fe_index, const unsigned int local_index) const { - Assert (fe_index == DoFHandler::default_fe_index, + Assert ((fe_index == DoFHandler::default_fe_index), ExcMessage ("Only the default FE index is allowed for non-hp DoFHandler objects")); Assert (selected_fe != 0, ExcMessage ("No finite element collection is associated with " @@ -1386,16 +1294,16 @@ get_vertex_dof_index (const unsigned int vertex_index, -template +template inline void -DoFHandler:: +DoFHandler:: set_vertex_dof_index (const unsigned int vertex_index, const unsigned int fe_index, const unsigned int local_index, const unsigned int global_index) { - Assert (fe_index == DoFHandler::default_fe_index, + Assert ((fe_index == DoFHandler::default_fe_index), ExcMessage ("Only the default FE index is allowed for non-hp DoFHandler objects")); Assert (selected_fe != 0, ExcMessage ("No finite element collection is associated with " @@ -1409,11 +1317,11 @@ set_vertex_dof_index (const unsigned int vertex_index, -template +template template inline bool -DoFHandler::fe_index_is_active (const unsigned int, +DoFHandler::fe_index_is_active (const unsigned int, const unsigned int, const unsigned int fe_index) const { @@ -1422,11 +1330,11 @@ DoFHandler::fe_index_is_active (const unsigned int, -template +template template inline unsigned int -DoFHandler::n_active_fe_indices (const unsigned int obj_level, +DoFHandler::n_active_fe_indices (const unsigned int obj_level, const unsigned int obj_index) const { // check that the object we look at is in @@ -1462,11 +1370,11 @@ DoFHandler::n_active_fe_indices (const unsigned int obj_level, -template +template template inline unsigned int -DoFHandler::nth_active_fe_index (const unsigned int obj_level, +DoFHandler::nth_active_fe_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int n) const { diff --git a/deal.II/deal.II/include/dofs/dof_iterator_selector.h b/deal.II/deal.II/include/dofs/dof_iterator_selector.h index b764033bae..807b19bdad 100644 --- a/deal.II/deal.II/include/dofs/dof_iterator_selector.h +++ b/deal.II/deal.II/include/dofs/dof_iterator_selector.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -13,165 +13,97 @@ #ifndef __deal2__dof_iterators_h #define __deal2__dof_iterators_h -//TODO: Rename this file to dof_iterators.h -//TODO: Can we avoid the specializations? - -// To consider for the second TODO: -// 1. Define accessors for unreasonable dimensions -// 2. Can cell_iterator and quad_iterator in 2D be different (or hex -// in 3D)? - #include + DEAL_II_NAMESPACE_OPEN +template class InvalidAccessor; + template class DoFAccessor; template class DoFCellAccessor; -template class DoFObjectAccessor; -template class DoFObjectAccessor<0,DH>; -template class DoFObjectAccessor<1,DH>; -template class DoFObjectAccessor<2,DH>; -template class DoFObjectAccessor<3,DH>; - -template class FiniteElement; -template class TriaRawIterator; -template class TriaIterator; -template class TriaActiveIterator; -template class Triangulation; -template class DoFHandler; +template class FiniteElement; +template class TriaRawIterator; +template class TriaIterator; +template class TriaActiveIterator; +template class Triangulation; +template class DoFHandler; namespace internal { namespace DoFHandler { -/** - * A pseudo class defining the iterator types used by DoFHandler and - * hp::DoFHandler. The typedefs in this class are synonymous with - * those in the handler classes, such that this class can be used in - * order to avoid including the complete header files of the handler. - * - * The class template is not used, but only its specializations with - * respect to the dimension. These differ from the template in two - * points: first, the iterator to objects of the same dimension as the - * handler is aliased to the #cell_iterator. Second, #face_iterator is - * aliased to objects of codimension one. - * - * - * @warning This raw iterators are not normally used in application - * programs. - * - * @ingroup dofs - * @ingroup Accessors - * @author W. Bangerth, G. Kanschat, O. Kayser-Herold, 1998, 2003, 2006 - */ template - struct Iterators - { - /// The dof handler class. - typedef DH DoFHandler_type; - /// The topological dimension of the dof handler. - static const unsigned int dim = DH::dimension; - - /// Iterator for raw lines. - typedef TriaRawIterator > raw_line_iterator; - /// Iterator for usual lines. - typedef TriaIterator > line_iterator; - /// Iterator for active lines. - typedef TriaActiveIterator > active_line_iterator; - /// Iterator for raw quadrilaterals - typedef TriaRawIterator > raw_quad_iterator; - /// Iterator for quadrilaterals - typedef TriaIterator > quad_iterator; - /// Iterator for active quadrilaterals - typedef TriaActiveIterator > active_quad_iterator; - /// Iterator for raw hexahedra - typedef TriaRawIterator > raw_hex_iterator; - /// Iterator for hexahedra - typedef TriaIterator > hex_iterator; - /// Iterator for active hexahedra - typedef TriaActiveIterator > active_hex_iterator; - /// Iterator for raw cells - typedef TriaRawIterator > raw_cell_iterator; - /// Iterator for cells - typedef TriaIterator > cell_iterator; - /// Iterator for active cells - typedef TriaActiveIterator > active_cell_iterator; - /// Iterator for raw faces - typedef TriaRawIterator > raw_face_iterator; - /// Iterator for faces - typedef TriaIterator > face_iterator; - /// Iterator for active faces - typedef TriaActiveIterator > active_face_iterator; - }; - + struct Iterators; /** - * Define some types for the DoF handling in one dimension. + * Define some types for DoF handling in one dimension. * * The types have the same meaning as those declared in - * internal::Triangulation::Iterators<2>, only the treatment of templates is a - * little more complicated since the @p DoFAccessor classes want a - * template template argument. + * internal::Triangulation::Iterators<1,spacedim>, only the treatment of + * templates is a little more complicated. See the @ref Iterators module + * for more information. * - * @author Wolfgang Bangerth, Oliver Kayser-Herold, 1998, 2003 + * @author Wolfgang Bangerth, Oliver Kayser-Herold, 1998, 2003, 2008 */ - template