]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Simplify get_dof_indices() and get_dof_indices() 10268/head
authorPeter Munch <peterrmuench@gmail.com>
Tue, 19 May 2020 10:35:51 +0000 (12:35 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Thu, 21 May 2020 20:01:48 +0000 (22:01 +0200)
include/deal.II/dofs/dof_accessor.h
include/deal.II/dofs/dof_accessor.templates.h
source/dofs/dof_accessor.cc
source/dofs/dof_handler_policy.cc

index 9f76c29489878af062a9db2282ccf8ed4276e9fb..c5ed6f50aa6886c02d099650334a1d5e5436ec13 100644 (file)
@@ -1288,6 +1288,16 @@ public:
   template <typename OtherAccessor>
   DoFInvalidAccessor(const OtherAccessor &);
 
+  /**
+   * Return the index of the <i>i</i>th 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<dim, spacedim>::default_fe_index) const;
+
   /**
    * Set the index of the <i>i</i>th degree of freedom of this object to @p
    * index. Since the current object doesn't point to anything useful, like
index af5eb4ce71c0e9e40f750aa948db237d6ddc68cb..4e1e865e3c4e8a14eb12eabb3f9c4ec14e9faeb1 100644 (file)
@@ -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 <typename DoFHandlerType,
+                bool level_dof_access,
+                int  structdim,
+                typename DoFIndicesType,
+                typename DoFOperation>
+      static void
+      process_dof_indices(
+        const dealii::DoFAccessor<structdim, DoFHandlerType, level_dof_access>
+          &                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<structdim>::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<structdim>::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<structdim>::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 <typename DoFHandlerType, bool level_dof_access, int structdim>
+      struct DoFIndexGetter
+      {
+        /**
+         * Return vertex DoF index.
+         */
+        static DEAL_II_ALWAYS_INLINE void
+        process_vertex_dof(
+          const dealii::DoFAccessor<structdim, DoFHandlerType, level_dof_access>
+            &                      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 <int structdim_>
+        static DEAL_II_ALWAYS_INLINE void
+        process_dof(const dealii::DoFAccessor<structdim_,
+                                              DoFHandlerType,
+                                              level_dof_access> &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 <int structdim_>
+        static DEAL_II_ALWAYS_INLINE void
+        process_dof(
+          const dealii::DoFInvalidAccessor<structdim_,
+                                           DoFHandlerType::dimension,
+                                           DoFHandlerType::space_dimension> &,
+          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 <typename DoFHandlerType, bool level_dof_access, int structdim>
+      struct DoFIndexSetter
+      {
+        /**
+         * Set vertex DoF index.
+         */
+        static DEAL_II_ALWAYS_INLINE void
+        process_vertex_dof(
+          const dealii::DoFAccessor<structdim, DoFHandlerType, level_dof_access>
+            &                            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 <int structdim_>
+        static DEAL_II_ALWAYS_INLINE void
+        process_dof(const dealii::DoFAccessor<structdim_,
+                                              DoFHandlerType,
+                                              level_dof_access> &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 <int structdim_>
+        static DEAL_II_ALWAYS_INLINE void
+        process_dof(
+          const dealii::DoFInvalidAccessor<structdim_,
+                                           DoFHandlerType::dimension,
+                                           DoFHandlerType::space_dimension> &,
+          const unsigned int,
+          const types::global_dof_index &,
+          const unsigned int)
+        {
+          Assert(false, ExcInternalError());
+        }
+      };
+
+
+
+      template <typename DoFHandlerType, bool level_dof_access, int structdim>
+      static void
+      get_dof_indices(
+        const dealii::DoFAccessor<structdim, DoFHandlerType, level_dof_access>
+          &                                   accessor,
+        std::vector<types::global_dof_index> &dof_indices,
+        const unsigned int                    fe_index)
+      {
+        process_dof_indices(
+          accessor,
+          dof_indices,
+          fe_index,
+          DoFIndexGetter<DoFHandlerType, level_dof_access, structdim>());
+      }
+
+
+
+      template <typename DoFHandlerType, bool level_dof_access, int structdim>
+      static void
+      set_dof_indices(
+        const dealii::DoFAccessor<structdim, DoFHandlerType, level_dof_access>
+          &                                         accessor,
+        const std::vector<types::global_dof_index> &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<DoFHandlerType, level_dof_access, structdim>());
+      }
     };
   } // namespace DoFAccessorImplementation
 } // namespace internal
@@ -1802,136 +2048,6 @@ namespace internal
 {
   namespace DoFAccessorImplementation
   {
-    template <typename DoFHandlerType, bool level_dof_access>
-    void
-    get_dof_indices(
-      const dealii::DoFAccessor<1, DoFHandlerType, level_dof_access> &accessor,
-      std::vector<types::global_dof_index> &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<types::global_dof_index>::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 <typename DoFHandlerType, bool level_dof_access>
-    void
-    get_dof_indices(
-      const dealii::DoFAccessor<2, DoFHandlerType, level_dof_access> &accessor,
-      std::vector<types::global_dof_index> &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<types::global_dof_index>::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 <typename DoFHandlerType, bool level_dof_access>
-    void
-    get_dof_indices(
-      const dealii::DoFAccessor<3, DoFHandlerType, level_dof_access> &accessor,
-      std::vector<types::global_dof_index> &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<types::global_dof_index>::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 <typename DoFHandlerType, bool level_dof_access>
     void
     get_mg_dof_indices(
@@ -2108,9 +2224,8 @@ DoFAccessor<structdim, DoFHandlerType, level_dof_access>::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 <typename DoFHandlerType, bool level_dof_access>
-      static void
-      set_dof_indices(
-        const DoFCellAccessor<DoFHandlerType, level_dof_access> &accessor,
-        const std::vector<types::global_dof_index> &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<dim>::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<dim>::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<dim>::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.
index 008758fadf8b47575fa3fdedb2a2016369c919e0..2b01478fdde9e269ca98772d25595fe88762825c 100644 (file)
@@ -75,6 +75,18 @@ DoFInvalidAccessor<structdim, dim, spacedim>::DoFInvalidAccessor(
 
 
 
+template <int structdim, int dim, int spacedim>
+types::global_dof_index
+DoFInvalidAccessor<structdim, dim, spacedim>::dof_index(
+  const unsigned int,
+  const unsigned int) const
+{
+  Assert(false, ExcInternalError());
+  return 0;
+}
+
+
+
 template <int structdim, int dim, int spacedim>
 void
 DoFInvalidAccessor<structdim, dim, spacedim>::set_dof_index(
@@ -118,8 +130,11 @@ DoFCellAccessor<DoFHandlerType, lda>::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<DoFHandlerType, lda, DoFHandlerType::dimension>(
+      *this, local_dof_indices, this->active_fe_index());
 }
 
 
index 9e40f700134a3e7ba0c76a5f27eb82cf0752e324..bca98b1ffd56597264a02f3137a2958a30e3108c 100644 (file)
@@ -1671,8 +1671,11 @@ namespace internal
                   std::vector<types::global_dof_index> 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)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.