From: Peter Munch Date: Tue, 19 May 2020 10:35:51 +0000 (+0200) Subject: Simplify get_dof_indices() and get_dof_indices() X-Git-Tag: v9.3.0-rc1~1559^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dd0f1d5403924e0b72984462f7d2cbb18d40647b;p=dealii.git Simplify get_dof_indices() and get_dof_indices() --- diff --git a/include/deal.II/dofs/dof_accessor.h b/include/deal.II/dofs/dof_accessor.h index 9f76c29489..c5ed6f50aa 100644 --- a/include/deal.II/dofs/dof_accessor.h +++ b/include/deal.II/dofs/dof_accessor.h @@ -1288,6 +1288,16 @@ public: template DoFInvalidAccessor(const OtherAccessor &); + /** + * Return the index of the ith degree of freedom of this object to + * @p index. Since the current object doesn't point to anything useful, like + * all other functions in this class this function only throws an exception. + */ + types::global_dof_index + dof_index(const unsigned int i, + const unsigned int fe_index = + DoFHandler::default_fe_index) const; + /** * Set the index of the ith degree of freedom of this object to @p * index. Since the current object doesn't point to anything useful, like diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index af5eb4ce71..4e1e865e3c 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -1588,6 +1588,252 @@ namespace internal *local_values_begin = read_write_vector[sorted_indices_pos[i]]; } #endif + + + + /** + * Loop over all degrees of freedom of the object described by the + * provided @p accessor and @p fe_index and perform the static functions + * provided by DoFOperation (set/get) on these. + */ + template + static void + process_dof_indices( + const dealii::DoFAccessor + & accessor, + DoFIndicesType & dof_indices, + const unsigned int fe_index, + const DoFOperation &) + { + 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; // + + const unsigned int inner_dofs = + structdim == 1 ? dofs_per_line : + (structdim == 2 ? dofs_per_quad : dofs_per_hex); + + unsigned int index = 0; + + // 1) VERTEX dofs + for (const unsigned int vertex : + GeometryInfo::vertex_indices()) + for (unsigned int d = 0; d < dofs_per_vertex; ++d, ++index) + DoFOperation::process_vertex_dof( + accessor, vertex, d, dof_indices[index], fe_index); + + // 2) 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/cell-local) ordering. + if (structdim == 2 || structdim == 3) + for (unsigned int line = 0; + line < GeometryInfo::lines_per_cell; + ++line) + for (unsigned int d = 0; d < dofs_per_line; ++d, ++index) + DoFOperation::process_dof( + *accessor.line(line), + accessor.get_fe(fe_index) + .adjust_line_dof_index_for_line_orientation( + d, accessor.line_orientation(line)), + dof_indices[index], + fe_index); + + // 3) 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 + if (structdim == 3) + for (unsigned int quad = 0; + quad < GeometryInfo::quads_per_cell; + ++quad) + for (unsigned int d = 0; d < dofs_per_quad; ++d, ++index) + DoFOperation::process_dof( + *accessor.quad(quad), + 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)), + dof_indices[index], + fe_index); + + // 4) INNER dofs + for (unsigned int d = 0; d < inner_dofs; ++d, ++index) + DoFOperation::process_dof(accessor, d, dof_indices[index], fe_index); + + AssertDimension(dof_indices.size(), index); + } + + + + /** + * An internal struct encapsulating the task of getting (vertex) + * DoF indices. + */ + template + struct DoFIndexGetter + { + /** + * Return vertex DoF index. + */ + static DEAL_II_ALWAYS_INLINE void + process_vertex_dof( + const dealii::DoFAccessor + & accessor, + const unsigned int vertex, + const unsigned int d, + types::global_dof_index &index_value, + const unsigned int fe_index) + { + index_value = accessor.vertex_dof_index(vertex, d, fe_index); + } + + /** + * Return DoF index for lines, quads, and inner degrees of freedom. + */ + template + static DEAL_II_ALWAYS_INLINE void + process_dof(const dealii::DoFAccessor &accessor, + const unsigned int d, + types::global_dof_index & index_value, + const unsigned int fe_index) + { + index_value = accessor.dof_index(d, fe_index); + } + + /** + * Fallback for DoFInvalidAccessor. + */ + template + static DEAL_II_ALWAYS_INLINE void + process_dof( + const dealii::DoFInvalidAccessor &, + const unsigned int, + types::global_dof_index &, + const unsigned int) + { + Assert(false, ExcInternalError()); + } + }; + + + + /** + * An internal struct encapsulating the task of setting (vertex) + * DoF indices. + */ + template + struct DoFIndexSetter + { + /** + * Set vertex DoF index. + */ + static DEAL_II_ALWAYS_INLINE void + process_vertex_dof( + const dealii::DoFAccessor + & accessor, + const unsigned int vertex, + const unsigned int d, + const types::global_dof_index &index_value, + const unsigned int fe_index) + { + accessor.set_vertex_dof_index(vertex, d, index_value, fe_index); + } + + /** + * Set DoF index for lines, quads, and inner degrees of freedom. + */ + template + static DEAL_II_ALWAYS_INLINE void + process_dof(const dealii::DoFAccessor &accessor, + const unsigned int d, + const types::global_dof_index & index_value, + const unsigned int fe_index) + { + accessor.set_dof_index(d, index_value, fe_index); + } + + /** + * Fallback for DoFInvalidAccessor. + */ + template + static DEAL_II_ALWAYS_INLINE void + process_dof( + const dealii::DoFInvalidAccessor &, + const unsigned int, + const types::global_dof_index &, + const unsigned int) + { + Assert(false, ExcInternalError()); + } + }; + + + + template + static void + get_dof_indices( + const dealii::DoFAccessor + & accessor, + std::vector &dof_indices, + const unsigned int fe_index) + { + process_dof_indices( + accessor, + dof_indices, + fe_index, + DoFIndexGetter()); + } + + + + template + static void + set_dof_indices( + const dealii::DoFAccessor + & accessor, + const std::vector &dof_indices, + const unsigned int fe_index) + { + // Note: this function is as general as `get_dof_indices()`. This + // assert is placed here since it is currently only used by the + // function DoFCellAccessor::set_dof_indices(), which is called by + // internal::DoFHandlerImplementation::Policy::Implementation::distribute_dofs(). + // In the case of new use cases, this assert can be removed. + Assert( + DoFHandlerType::dimension == structdim, + ExcMessage( + "This function is intended to be used for DoFCellAccessor, i.e., dimension == structdim.")); + + process_dof_indices( + accessor, + dof_indices, + fe_index, + DoFIndexSetter()); + } }; } // namespace DoFAccessorImplementation } // namespace internal @@ -1802,136 +2048,6 @@ namespace internal { namespace DoFAccessorImplementation { - template - void - get_dof_indices( - const dealii::DoFAccessor<1, DoFHandlerType, level_dof_access> &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 : GeometryInfo<1>::vertex_indices()) - for (unsigned int d = 0; d < dofs_per_vertex; ++d) - *next++ = accessor.vertex_dof_index(vertex, d, fe_index); - for (unsigned int d = 0; d < dofs_per_line; ++d) - *next++ = accessor.dof_index(d, fe_index); - } - - - - template - void - get_dof_indices( - const dealii::DoFAccessor<2, DoFHandlerType, level_dof_access> &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 : GeometryInfo<2>::vertex_indices()) - for (unsigned int d = 0; d < dofs_per_vertex; ++d) - *next++ = accessor.vertex_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 < GeometryInfo<2>::lines_per_cell; - ++line) - for (unsigned int d = 0; d < dofs_per_line; ++d) - *next++ = accessor.line(line)->dof_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 < dofs_per_quad; ++d) - *next++ = accessor.dof_index(d, fe_index); - } - - - - template - void - get_dof_indices( - const dealii::DoFAccessor<3, DoFHandlerType, level_dof_access> &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 : GeometryInfo<3>::vertex_indices()) - for (unsigned int d = 0; d < dofs_per_vertex; ++d) - *next++ = accessor.vertex_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 < GeometryInfo<3>::lines_per_cell; - ++line) - for (unsigned int d = 0; d < dofs_per_line; ++d) - *next++ = accessor.line(line)->dof_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 < GeometryInfo<3>::faces_per_cell; - ++quad) - for (unsigned int d = 0; d < dofs_per_quad; ++d) - *next++ = accessor.quad(quad)->dof_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 < dofs_per_hex; ++d) - *next++ = accessor.dof_index(d, fe_index); - } - - - template void get_mg_dof_indices( @@ -2108,9 +2224,8 @@ DoFAccessor::get_dof_indices( ExcInternalError()); // now do the actual work - dealii::internal::DoFAccessorImplementation::get_dof_indices(*this, - dof_indices, - fe_index); + dealii::internal::DoFAccessorImplementation::Implementation::get_dof_indices( + *this, dof_indices, fe_index); } @@ -2715,88 +2830,6 @@ namespace internal - /** - * Implement setting dof indices on a cell. - */ - template - static void - set_dof_indices( - const DoFCellAccessor &accessor, - const std::vector &local_dof_indices) - { - const unsigned int dim = DoFHandlerType::dimension; - - Assert(accessor.has_children() == false, ExcInternalError()); - - 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; - - Assert(local_dof_indices.size() == accessor.get_fe().dofs_per_cell, - ExcInternalError()); - - unsigned int index = 0; - - for (const unsigned int vertex : GeometryInfo::vertex_indices()) - for (unsigned int d = 0; d < dofs_per_vertex; ++d, ++index) - accessor.set_vertex_dof_index(vertex, - d, - local_dof_indices[index], - accessor.active_fe_index()); - // now copy dof numbers into the line. for lines in 3d 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. - // - // of course, if dim<3, then there is nothing to adjust - for (unsigned int line = 0; line < GeometryInfo::lines_per_cell; - ++line) - for (unsigned int d = 0; d < dofs_per_line; ++d, ++index) - accessor.line(line)->set_dof_index( - dim < 3 ? - d : - accessor.get_fe().adjust_line_dof_index_for_line_orientation( - d, accessor.line_orientation(line)), - local_dof_indices[index], - accessor.active_fe_index()); - // now copy dof numbers into the face. for faces in 3d 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 - // - // again, if dim<3, then there is nothing to adjust - for (unsigned int quad = 0; quad < GeometryInfo::quads_per_cell; - ++quad) - for (unsigned int d = 0; d < dofs_per_quad; ++d, ++index) - accessor.quad(quad)->set_dof_index( - dim < 3 ? - d : - accessor.get_fe().adjust_quad_dof_index_for_face_orientation( - d, - accessor.face_orientation(quad), - accessor.face_flip(quad), - accessor.face_rotation(quad)), - local_dof_indices[index], - accessor.active_fe_index()); - for (unsigned int d = 0; d < dofs_per_hex; ++d, ++index) - accessor.set_dof_index(d, - local_dof_indices[index], - accessor.active_fe_index()); - - Assert(index == accessor.get_fe().dofs_per_cell, ExcInternalError()); - } - - - /** * Do what the active_fe_index function in the parent class is supposed to * do. diff --git a/source/dofs/dof_accessor.cc b/source/dofs/dof_accessor.cc index 008758fadf..2b01478fdd 100644 --- a/source/dofs/dof_accessor.cc +++ b/source/dofs/dof_accessor.cc @@ -75,6 +75,18 @@ DoFInvalidAccessor::DoFInvalidAccessor( +template +types::global_dof_index +DoFInvalidAccessor::dof_index( + const unsigned int, + const unsigned int) const +{ + Assert(false, ExcInternalError()); + return 0; +} + + + template void DoFInvalidAccessor::set_dof_index( @@ -118,8 +130,11 @@ DoFCellAccessor::set_dof_indices( Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject()); - internal::DoFCellAccessorImplementation::Implementation::set_dof_indices( - *this, local_dof_indices); + AssertDimension(local_dof_indices.size(), this->get_fe().dofs_per_cell); + + internal::DoFAccessorImplementation::Implementation:: + template set_dof_indices( + *this, local_dof_indices, this->active_fe_index()); } diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc index 9e40f70013..bca98b1ffd 100644 --- a/source/dofs/dof_handler_policy.cc +++ b/source/dofs/dof_handler_policy.cc @@ -1671,8 +1671,11 @@ namespace internal std::vector dof_indices( cell->get_fe().dofs_per_cell); - internal::DoFAccessorImplementation::get_dof_indices( - *cell, dof_indices, cell->active_fe_index()); + // circumvent cache + internal::DoFAccessorImplementation::Implementation:: + get_dof_indices(*cell, + dof_indices, + cell->active_fe_index()); for (auto &dof_index : dof_indices) if (dof_index == numbers::invalid_dof_index)