From: Wolfgang Bangerth Date: Thu, 30 Mar 2006 16:12:01 +0000 (+0000) Subject: Implement a cache for the results of get_dof_indices X-Git-Tag: v8.0.0~11942 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=408a0bc5474945b7f332a6f3ba54568333149aa0;p=dealii.git Implement a cache for the results of get_dof_indices git-svn-id: https://svn.dealii.org/trunk@12778 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_accessor.h b/deal.II/deal.II/include/dofs/dof_accessor.h index 9302b344e1..56915ac1b4 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.h @@ -329,6 +329,22 @@ class DoFObjectAccessor : public DoFAccessor, * 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. */ void get_dof_indices (std::vector &dof_indices) const; @@ -649,6 +665,22 @@ class DoFObjectAccessor<1, DH> : * * It is assumed that the vector already * has the right size beforehand. + * + * 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. */ void get_dof_indices (std::vector &dof_indices) const; @@ -896,6 +928,22 @@ class DoFObjectAccessor<2, DH> : * * It is assumed that the vector already * has the right size beforehand. + * + * 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. */ void get_dof_indices (std::vector &dof_indices) const; @@ -1151,6 +1199,22 @@ class DoFObjectAccessor<3, DH> : * * It is assumed that the vector already * has the right size beforehand. + * + * 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. */ void get_dof_indices (std::vector &dof_indices) const; @@ -1583,7 +1647,69 @@ class DoFCellAccessor : public DoFObjectAccessor template void set_dof_values_by_interpolation (const Vector &local_values, OutputVector &values) const; + + /** + * Return the indices of the dofs of this + * quad 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. + * + * It is assumed that the vector already + * has the right size beforehand. + * + * This function reimplements the + * same function in the base + * class. The functions in the + * base classes are available for + * all geometric objects, + * i.e. even in 3d they can be + * used to access the dof indices + * of edges, for example. On the + * other hand, the most common + * case is clearly the use on + * cells, which is why we cache + * the array for each cell, but + * not edge. To retrieve the + * cached values, rather than + * collect the necessary + * information every time, this + * function overwrites the one in + * the base class. + * + * 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. + */ + void get_dof_indices (std::vector &dof_indices) const; + private: + /** + * Update the cache in which we + * store the dof indices of this + * cell. + */ + void update_cell_dof_indices_cache () const; + + /** + * Make the DoFHandler class a + * friend so that it can call the + * update_cell_dof_indices_cache() + * function + */ + template friend class DoFHandler; }; @@ -1677,6 +1803,28 @@ void DoFObjectAccessor<3,hp::DoFHandler<3> >::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; + + + // include more templates #include "dof_accessor.templates.h" 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 b114d7833d..5c5d63f909 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.templates.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.templates.h @@ -134,7 +134,7 @@ unsigned int DoFObjectAccessor<1,DH>::dof_index (const unsigned int i) const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // since the exception classes are // from a template dependent base @@ -166,7 +166,7 @@ DoFObjectAccessor<1,DH>::vertex_dof_index (const unsigned int vertex, const unsigned int i) const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // since the exception classes are // from a template dependent base @@ -198,7 +198,7 @@ void DoFObjectAccessor<1,DH>::get_dof_indices (std::vector &dof_indices) const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // since the exception classes are // from a template dependent base @@ -248,7 +248,7 @@ TriaIterator > DoFObjectAccessor<1,DH>::child (const unsigned int i) const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); TriaIterator > q (this->tria, this->present_level+1, @@ -451,7 +451,7 @@ DoFObjectAccessor<2,DH>::dof_index (const unsigned int i) const Assert (iget_fe().dofs_per_quad, ExcIndexRange (i, 0, this->get_fe().dofs_per_quad)); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); return this->dof_handler->levels[this->present_level] ->quad_dofs[this->present_index*this->get_fe().dofs_per_quad+i]; @@ -474,7 +474,7 @@ DoFObjectAccessor<2,DH>::vertex_dof_index (const unsigned int vertex, Assert (iget_fe().dofs_per_vertex, ExcIndexRange (i, 0, this->get_fe().dofs_per_vertex)); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); const unsigned int dof_number = (this->vertex_index(vertex) * this->get_fe().dofs_per_vertex + @@ -483,7 +483,6 @@ DoFObjectAccessor<2,DH>::vertex_dof_index (const unsigned int vertex, } - template inline void @@ -514,7 +513,7 @@ DoFObjectAccessor<2,DH>::get_dof_indices (std::vector &dof_indices typename BaseClass::ExcNotActive()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); const unsigned int dofs_per_vertex = this->get_fe().dofs_per_vertex, dofs_per_line = this->get_fe().dofs_per_line, @@ -538,7 +537,7 @@ DoFObjectAccessor<2,DH>::line (const unsigned int i) const { Assert (i<4, ExcIndexRange (i, 0, 4)); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); return TriaIterator > ( @@ -557,7 +556,7 @@ TriaIterator > DoFObjectAccessor<2,DH>::child (const unsigned int i) const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); TriaIterator > q (this->tria, this->present_level+1, @@ -635,7 +634,7 @@ distribute_local_to_global (const Vector &local_source, typename BaseClass::ExcVectorDoesNotMatch()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); const unsigned int n_dofs = local_source.size(); @@ -686,7 +685,7 @@ distribute_local_to_global (const FullMatrix &local_source, Assert (global_destination.m() == global_destination.n(), typename BaseClass::ExcMatrixDoesNotMatch()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); const unsigned int n_dofs = local_source.m(); @@ -772,7 +771,7 @@ DoFObjectAccessor<3,DH>::dof_index (const unsigned int i) const Assert (iget_fe().dofs_per_hex, ExcIndexRange (i, 0, this->get_fe().dofs_per_hex)); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); return this->dof_handler->levels[this->present_level] ->hex_dofs[this->present_index*this->get_fe().dofs_per_hex+i]; @@ -796,7 +795,7 @@ DoFObjectAccessor<3,DH>::vertex_dof_index (const unsigned int vertex, Assert (iget_fe().dofs_per_vertex, ExcIndexRange (i, 0, this->get_fe().dofs_per_vertex)); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); const unsigned int dof_number = (this->vertex_index(vertex) * this->get_fe().dofs_per_vertex + @@ -822,7 +821,7 @@ DoFObjectAccessor<3,DH>::get_dof_indices (std::vector &dof_indices this->get_fe().dofs_per_hex), typename BaseClass::ExcVectorDoesNotMatch()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // this function really only makes // sense on non-active objects if @@ -837,7 +836,7 @@ DoFObjectAccessor<3,DH>::get_dof_indices (std::vector &dof_indices 8*this->get_fe().dofs_per_vertex), typename BaseClass::ExcNotActive()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); const unsigned int dofs_per_vertex = this->get_fe().dofs_per_vertex, dofs_per_line = this->get_fe().dofs_per_line, @@ -969,7 +968,7 @@ distribute_local_to_global (const Vector &local_source, Assert (this->dof_handler->n_dofs() == global_destination.size(), typename BaseClass::ExcVectorDoesNotMatch()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); const unsigned int n_dofs = local_source.size(); @@ -1021,7 +1020,7 @@ distribute_local_to_global (const FullMatrix &local_source, Assert (global_destination.m() == global_destination.n(), typename BaseClass::ExcMatrixDoesNotMatch()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); const unsigned int n_dofs = local_source.m(); @@ -1090,7 +1089,7 @@ DoFObjectAccessor<1,hp::DoFHandler<1> >::dof_index (const unsigned int i) const Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, BaseClass::ExcInvalidObject()); @@ -1127,7 +1126,7 @@ DoFObjectAccessor<1,hp::DoFHandler<2> >::dof_index (const unsigned int i) const Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, BaseClass::ExcInvalidObject()); @@ -1164,7 +1163,7 @@ DoFObjectAccessor<1,hp::DoFHandler<3> >::dof_index (const unsigned int i) const Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, BaseClass::ExcInvalidObject()); @@ -1195,7 +1194,7 @@ DoFObjectAccessor<1,hp::DoFHandler<1> >::set_dof_index (const unsigned int i, Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, @@ -1226,7 +1225,7 @@ DoFObjectAccessor<1, hp::DoFHandler<2> >::set_dof_index (const unsigned int i, Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, @@ -1261,7 +1260,7 @@ DoFObjectAccessor<1, hp::DoFHandler<3> >::set_dof_index (const unsigned int i, Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, @@ -1319,7 +1318,7 @@ unsigned int DoFObjectAccessor<1,hp::DoFHandler<1> >::active_fe_index () const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1336,7 +1335,7 @@ unsigned int DoFObjectAccessor<1,hp::DoFHandler<2> >::active_fe_index () const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1353,7 +1352,7 @@ unsigned int DoFObjectAccessor<1,hp::DoFHandler<3> >::active_fe_index () const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1372,7 +1371,7 @@ DoFObjectAccessor<1,hp::DoFHandler<1> >::set_active_fe_index (const unsigned int Assert (this->dof_handler != 0, ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1392,7 +1391,7 @@ DoFObjectAccessor<1,hp::DoFHandler<2> >::set_active_fe_index (const unsigned int { Assert (this->dof_handler != 0, ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1413,7 +1412,7 @@ DoFObjectAccessor<1,hp::DoFHandler<3> >::set_active_fe_index (const unsigned int { Assert (this->dof_handler != 0, ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1437,7 +1436,7 @@ unsigned int DoFObjectAccessor<2,hp::DoFHandler<2> >::dof_index (const unsigned Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, @@ -1467,7 +1466,7 @@ unsigned int DoFObjectAccessor<2,hp::DoFHandler<3> >::dof_index (const unsigned Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, @@ -1499,7 +1498,7 @@ DoFObjectAccessor<2, hp::DoFHandler<2> >::set_dof_index (const unsigned int i, Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, @@ -1531,7 +1530,7 @@ DoFObjectAccessor<2, hp::DoFHandler<3> >::set_dof_index (const unsigned int i, Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, @@ -1578,7 +1577,7 @@ unsigned int DoFObjectAccessor<2,hp::DoFHandler<2> >::active_fe_index () const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1595,7 +1594,7 @@ unsigned int DoFObjectAccessor<2,hp::DoFHandler<3> >::active_fe_index () const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1612,9 +1611,9 @@ DoFObjectAccessor<2,hp::DoFHandler<2> >::set_active_fe_index (const unsigned int { Assert (this->dof_handler != 0, ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1635,7 +1634,7 @@ DoFObjectAccessor<2,hp::DoFHandler<3> >::set_active_fe_index (const unsigned int { Assert (this->dof_handler != 0, ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1660,7 +1659,7 @@ DoFObjectAccessor<3,hp::DoFHandler<3> >::dof_index (const unsigned int i) const Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, @@ -1691,7 +1690,7 @@ DoFObjectAccessor<3, hp::DoFHandler<3> >::set_dof_index (const unsigned int i, Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); // make sure a FE has been selected // and enough room was reserved Assert (&this->get_fe() != 0, @@ -1728,7 +1727,7 @@ unsigned int DoFObjectAccessor<3,hp::DoFHandler<3> >::active_fe_index () const { Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1746,7 +1745,7 @@ DoFObjectAccessor<3,hp::DoFHandler<3> >::set_active_fe_index (const unsigned int { Assert (this->dof_handler != 0, ExcInvalidObject()); Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), - ExcMessage ("hp::DoFHandler not initialized")); + ExcMessage ("DoFHandler not initialized")); Assert (static_cast::size_type>(this->present_index) < this->dof_handler->levels[this->present_level]->active_fe_indices.size (), ExcIndexRange (this->present_index, 0, @@ -1840,4 +1839,119 @@ DoFCellAccessor >::face (const unsigned int i) const } + + +template <> +inline +void +DoFCellAccessor >:: +get_dof_indices (std::vector &dof_indices) const +{ + Assert (dof_indices.size() == this->get_fe().dofs_per_cell, + 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->present_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); +} + + + +template <> +inline +void +DoFCellAccessor >:: +get_dof_indices (std::vector &dof_indices) const +{ + Assert (dof_indices.size() == this->get_fe().dofs_per_cell, + 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->present_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); +} + + + +template <> +inline +void +DoFCellAccessor >:: +get_dof_indices (std::vector &dof_indices) const +{ + Assert (dof_indices.size() == this->get_fe().dofs_per_cell, + 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<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->present_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); +} + + #endif diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index 2a4af7d8b7..9d5823292f 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -1055,6 +1055,11 @@ class DoFHandler : public Subscriptor * Make accessor objects friends. */ template friend class DoFObjectAccessor; + + /* + * Same with cell accessor + */ + template friend class DoFCellAccessor; }; diff --git a/deal.II/deal.II/include/dofs/dof_levels.h b/deal.II/deal.II/include/dofs/dof_levels.h index bbfb1f994e..a361c9a0fd 100644 --- a/deal.II/deal.II/include/dofs/dof_levels.h +++ b/deal.II/deal.II/include/dofs/dof_levels.h @@ -46,17 +46,24 @@ namespace internal * * The storage format of the degrees of freedom indices (short: DoF * indices) is somewhat like a mirror of the data structures of the - * triangulation classes. There is a hierarchy of - * DoFLevel classes for the different dimensions which - * have objects named @p line_dofs, @p quad_dofs and so on, in which - * the indices of DoFs located on lines and quads, respectively, are - * stored. The indices are stored levelwise. The layout in these - * arrays is as follows: if for a selected finite element (use - * DoFHandler::distribute_dofs() to select a finite element) - * the number of DoFs on each line (without those in the vertices) is - * @p N, then the length of the @p line_dofs array is @p N times - * the number of lines on this level. The DoF indices for the @p ith - * line are at the positions N*i...(N+1)*i-1. + * triangulation classes. There is a hierarchy of DoFLevel + * classes for the different dimensions which have objects named @p + * line_dofs, @p quad_dofs and so on, in which the indices of DoFs + * located on lines and quads, respectively, are stored. The indices + * are stored levelwise. The layout in these arrays is as follows: if + * for a selected finite element (use DoFHandler::distribute_dofs() to + * select a finite element) the number of DoFs on each line (without + * those in the vertices) is @p N, then the length of the @p line_dofs + * array is @p N times the number of lines on this level. The DoF + * indices for the @p ith line are at the positions + * N*i...(N+1)*i-1. It is true that this scheme wastes some + * memory since we actually only use the memory so allocated for + * active lines, whereas the index entries for line i remain unused if + * this line has children. On the other hand, it is easy to see that + * at least for dim>1 does not constitute a too big fraction of memory + * entries, and we opt for higher memory consumption in return for + * less index computations that usually slow down modern processors + * due to more complex memory access patterns. * * Since vertices are not associated with a particular level, the * indices associated with vertices are not stored in the DoFLevel @@ -78,6 +85,43 @@ namespace internal }; + +/** + * Common base class for all the DoFLevel classes. We here store + * information that is associated with (logical) cells, rather than + * concrete objects such as lines, quads, or hexes. + * + * At present, all we store are cached values for the DoF indices on + * each cell, since this is a frequently requested operation. The + * values are set by DoFCellAccessor::update_cell_dof_indices_cache + * and are used by DoFCellAccessor::get_dof_indices. Note that for + * objects that are not cell, such as lines in 2d and 3d, or quads in + * 3d, no caching is done since accessing the degrees of freedom on + * these objects is a comparatively rare operation. + */ + template <> + class DoFLevel<0> + { + public: + /** + * Cache for the DoF indices + * on cells. The size of this + * array equals the number of + * cells on a given level + * times + * selected_fe.dofs_per_cell. + */ + std::vector cell_dof_indices_cache; + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. + */ + unsigned int memory_consumption () const; + }; + + /** * Store the indices of the degrees of freedom which are located on * the lines. See the general template DoFLevel for more information. @@ -86,7 +130,7 @@ namespace internal * @author Wolfgang Bangerth, 1998 */ template <> - class DoFLevel<1> + class DoFLevel<1> : public DoFLevel<0> { public: /** diff --git a/deal.II/deal.II/source/dofs/dof_accessor.cc b/deal.II/deal.II/source/dofs/dof_accessor.cc index a3f502d19b..8c57ae088c 100644 --- a/deal.II/deal.II/source/dofs/dof_accessor.cc +++ b/deal.II/deal.II/source/dofs/dof_accessor.cc @@ -432,6 +432,68 @@ DoFCellAccessor >::face (const unsigned int) const return TriaIterator<1, DoFObjectAccessor<0,DoFHandler<1> > >(); } + +template <> +void +DoFCellAccessor >::update_cell_dof_indices_cache () const +{ + Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), + ExcMessage ("DoFHandler not initialized")); + + Assert (this->dof_handler != 0, DoFAccessor >::ExcInvalidObject()); + Assert (&this->get_fe() != 0, DoFAccessor >::ExcInvalidObject()); + + // 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 (this->has_children() + && + (this->get_fe().dofs_per_cell != + this->get_fe().dofs_per_vertex * GeometryInfo<1>::vertices_per_cell)) + return; + + const unsigned int dofs_per_vertex = this->get_fe().dofs_per_vertex, + dofs_per_line = this->get_fe().dofs_per_line, + dofs_per_cell = this->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 (this->present_index * dofs_per_cell + dofs_per_cell + <= + this->dof_handler->levels[this->present_level] + ->cell_dof_indices_cache.size(), + ExcInternalError()); + + std::vector::iterator next + = (this->dof_handler->levels[this->present_level] + ->cell_dof_indices_cache.begin() + this->present_index * dofs_per_cell); + + for (unsigned int vertex=0; vertex<2; ++vertex) + for (unsigned int d=0; d +void +DoFCellAccessor >::update_cell_dof_indices_cache () const +{ + // not implemented, but should also + // not be called + Assert (false, ExcNotImplemented()); +} + + + #endif @@ -444,6 +506,72 @@ DoFCellAccessor >::face (const unsigned int i) const return this->line(i); } + + +template <> +void +DoFCellAccessor >::update_cell_dof_indices_cache () const +{ + Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), + ExcMessage ("DoFHandler not initialized")); + + Assert (this->dof_handler != 0, DoFAccessor >::ExcInvalidObject()); + Assert (&this->get_fe() != 0, DoFAccessor >::ExcInvalidObject()); + + // 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 (this->has_children() + && + (this->get_fe().dofs_per_cell != + this->get_fe().dofs_per_vertex * GeometryInfo<2>::vertices_per_cell)) + return; + + const unsigned int dofs_per_vertex = this->get_fe().dofs_per_vertex, + dofs_per_line = this->get_fe().dofs_per_line, + dofs_per_quad = this->get_fe().dofs_per_quad, + dofs_per_cell = this->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 (this->present_index * dofs_per_cell + dofs_per_cell + <= + this->dof_handler->levels[this->present_level] + ->cell_dof_indices_cache.size(), + ExcInternalError()); + + std::vector::iterator next + = (this->dof_handler->levels[this->present_level] + ->cell_dof_indices_cache.begin() + this->present_index * dofs_per_cell); + + for (unsigned int vertex=0; vertex<4; ++vertex) + for (unsigned int d=0; dline(line)->dof_index(d); + for (unsigned int d=0; d +void +DoFCellAccessor >::update_cell_dof_indices_cache () const +{ + // not implemented, but should also + // not be called + Assert (false, ExcNotImplemented()); +} + + #endif @@ -456,6 +584,75 @@ DoFCellAccessor >::face (const unsigned int i) const return this->quad(i); } + + +template <> +void +DoFCellAccessor >::update_cell_dof_indices_cache () const +{ + Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), + ExcMessage ("DoFHandler not initialized")); + + Assert (this->dof_handler != 0, DoFAccessor >::ExcInvalidObject()); + Assert (&this->get_fe() != 0, DoFAccessor >::ExcInvalidObject()); + + // 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 (this->has_children() + && + (this->get_fe().dofs_per_cell != + this->get_fe().dofs_per_vertex * GeometryInfo<3>::vertices_per_cell)) + return; + + const unsigned int dofs_per_vertex = this->get_fe().dofs_per_vertex, + dofs_per_line = this->get_fe().dofs_per_line, + dofs_per_quad = this->get_fe().dofs_per_quad, + dofs_per_hex = this->get_fe().dofs_per_hex, + dofs_per_cell = this->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 (this->present_index * dofs_per_cell + dofs_per_cell + <= + this->dof_handler->levels[this->present_level] + ->cell_dof_indices_cache.size(), + ExcInternalError()); + + std::vector::iterator next + = (this->dof_handler->levels[this->present_level] + ->cell_dof_indices_cache.begin() + this->present_index * dofs_per_cell); + + for (unsigned int vertex=0; vertex<8; ++vertex) + for (unsigned int d=0; dline(line)->dof_index(d); + for (unsigned int quad=0; quad<6; ++quad) + for (unsigned int d=0; dquad(quad)->dof_index(d); + for (unsigned int d=0; d +void +DoFCellAccessor >::update_cell_dof_indices_cache () const +{ + // not implemented, but should also + // not be called + Assert (false, ExcNotImplemented()); +} + + #endif diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index be621d115c..032b5e8afb 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -1457,13 +1457,20 @@ void DoFHandler::distribute_dofs (const FiniteElement &ff, used_dofs = next_free_dof; + // update the cache used + // for cell dof indices + for (cell_iterator cell = begin(); cell != end(); ++cell) + cell->update_cell_dof_indices_cache (); + // finally restore the user flags const_cast &>(*tria).load_user_flags(user_flags); } + template -void DoFHandler::clear () { +void DoFHandler::clear () +{ // release lock to old fe selected_fe = 0; @@ -1682,6 +1689,11 @@ void DoFHandler<1>::renumber_dofs (const std::vector &new_numbers) i!=levels[level]->line_dofs.end(); ++i) if (*i != invalid_dof_index) *i = new_numbers[*i]; + + // finally update the cache used + // for cell dof indices + for (cell_iterator cell = begin(); cell != end(); ++cell) + cell->update_cell_dof_indices_cache (); } #endif @@ -1735,7 +1747,12 @@ void DoFHandler<2>::renumber_dofs (const std::vector &new_numbers) i!=levels[level]->quad_dofs.end(); ++i) if (*i != invalid_dof_index) *i = new_numbers[*i]; - }; + } + + // finally update the cache used + // for cell dof indices + for (cell_iterator cell = begin(); cell != end(); ++cell) + cell->update_cell_dof_indices_cache (); } #endif @@ -1793,7 +1810,12 @@ void DoFHandler<3>::renumber_dofs (const std::vector &new_numbers) i!=levels[level]->hex_dofs.end(); ++i) if (*i != invalid_dof_index) *i = new_numbers[*i]; - }; + } + + // finally update the cache used + // for cell dof indices + for (cell_iterator cell = begin(); cell != end(); ++cell) + cell->update_cell_dof_indices_cache (); } #endif @@ -1967,7 +1989,8 @@ DoFHandler<3>::max_couplings_between_boundary_dofs () const #if deal_II_dimension == 1 template <> -void DoFHandler<1>::reserve_space () { +void DoFHandler<1>::reserve_space () +{ Assert (selected_fe != 0, ExcNoFESelected()); Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); @@ -1984,10 +2007,14 @@ void DoFHandler<1>::reserve_space () { { levels.push_back (new internal::DoFHandler::DoFLevel<1>); - levels.back()->line_dofs = std::vector(tria->n_raw_lines(i) * - selected_fe->dofs_per_line, - invalid_dof_index); - }; + levels.back()->line_dofs.resize (tria->n_raw_lines(i) * + selected_fe->dofs_per_line, + invalid_dof_index); + + levels.back()->cell_dof_indices_cache.resize (tria->n_raw_lines(i) * + selected_fe->dofs_per_cell, + invalid_dof_index); + } } @@ -1997,7 +2024,8 @@ void DoFHandler<1>::reserve_space () { #if deal_II_dimension == 2 template <> -void DoFHandler<2>::reserve_space () { +void DoFHandler<2>::reserve_space () +{ Assert (selected_fe != 0, ExcNoFESelected()); Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); @@ -2014,13 +2042,17 @@ void DoFHandler<2>::reserve_space () { { levels.push_back (new internal::DoFHandler::DoFLevel<2>); - levels.back()->line_dofs = std::vector (tria->n_raw_lines(i) * - selected_fe->dofs_per_line, - invalid_dof_index); - levels.back()->quad_dofs = std::vector (tria->n_raw_quads(i) * - selected_fe->dofs_per_quad, - invalid_dof_index); - }; + levels.back()->line_dofs.resize (tria->n_raw_lines(i) * + selected_fe->dofs_per_line, + invalid_dof_index); + levels.back()->quad_dofs.resize (tria->n_raw_quads(i) * + selected_fe->dofs_per_quad, + invalid_dof_index); + + levels.back()->cell_dof_indices_cache.resize (tria->n_raw_quads(i) * + selected_fe->dofs_per_cell, + invalid_dof_index); + } } #endif @@ -2029,7 +2061,8 @@ void DoFHandler<2>::reserve_space () { #if deal_II_dimension == 3 template <> -void DoFHandler<3>::reserve_space () { +void DoFHandler<3>::reserve_space () +{ Assert (selected_fe != 0, ExcNoFESelected()); Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); @@ -2055,7 +2088,11 @@ void DoFHandler<3>::reserve_space () { levels.back()->hex_dofs = std::vector (tria->n_raw_hexs(i) * selected_fe->dofs_per_hex, invalid_dof_index); - }; + + levels.back()->cell_dof_indices_cache.resize (tria->n_raw_hexs(i) * + selected_fe->dofs_per_cell, + invalid_dof_index); + } } #endif diff --git a/deal.II/deal.II/source/dofs/dof_levels.cc b/deal.II/deal.II/source/dofs/dof_levels.cc index 3cf22e4417..2f97a2661b 100644 --- a/deal.II/deal.II/source/dofs/dof_levels.cc +++ b/deal.II/deal.II/source/dofs/dof_levels.cc @@ -21,10 +21,18 @@ namespace internal namespace DoFHandler { + unsigned int + DoFLevel<0>::memory_consumption () const + { + return MemoryConsumption::memory_consumption (cell_dof_indices_cache); + } + + unsigned int DoFLevel<1>::memory_consumption () const { - return MemoryConsumption::memory_consumption (line_dofs); + return (DoFLevel<0>::memory_consumption () + + MemoryConsumption::memory_consumption (line_dofs)); } diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index 6cc91c7e76..d48dad6feb 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -552,6 +552,15 @@ inconvenience this causes.

deal.II

    +

  1. + New: Since calling cell->get_dof_indices is a + fairly frequent operation (called about 7 times per cell in + step-14, but much more often in more complicated programs), the + results of this function are now cached for faster access. +
    + (WB 2006/03/29) +

    +

  2. Fixed: An exception was generated when trying to interpolate boundary values for a primitive component of a composite finite