]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Compute multigrid dof handler also for distributed policy 3D
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 22 Sep 2015 16:15:54 +0000 (18:15 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 23 Sep 2015 15:02:22 +0000 (17:02 +0200)
include/deal.II/dofs/dof_accessor.templates.h
source/dofs/dof_handler_policy.cc

index d405829cb2e3b9ceb7c5f3221997e8ba29f2d054..f879e4e0f909d5c4c3bddb3ea61b7ea63bb81527 100644 (file)
@@ -1433,12 +1433,30 @@ namespace internal
 
       template<class DH, bool level_dof_access>
       static
-      void set_mg_dof_indices (const dealii::DoFAccessor<3, DH,level_dof_access> &,
-                               const int,
-                               const std::vector<types::global_dof_index> &,
-                               const unsigned int)
+      void set_mg_dof_indices (const dealii::DoFAccessor<3, DH,level_dof_access> &accessor,
+                               const int level,
+                               const std::vector<types::global_dof_index> &dof_indices,
+                               const unsigned int fe_index)
       {
-        AssertThrow (false, ExcNotImplemented ()); //TODO[TH]: implement
+        const FiniteElement<DH::dimension, DH::space_dimension> &fe = accessor.get_dof_handler ().get_fe ()[fe_index];
+        std::vector<types::global_dof_index>::const_iterator next = dof_indices.begin ();
+
+        for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex)
+          for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
+            accessor.set_mg_vertex_dof_index(level, vertex, dof, *next++, fe_index);
+
+        for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell; ++line)
+          for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
+            accessor.line(line)->set_mg_dof_index(level, dof, *next++);
+
+        for (unsigned int quad = 0; quad < GeometryInfo<3>::quads_per_cell; ++quad)
+          for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof)
+            accessor.quad(quad)->set_mg_dof_index(level, dof, *next++);
+
+        for (unsigned int dof = 0; dof < fe.dofs_per_hex; ++dof)
+          accessor.set_mg_dof_index(level, dof, *next++);
+
+        Assert (next == dof_indices.end (), ExcInternalError ());
       }
 
     };
index dd3eaff84e8eb3d650959627165dcae443c863d6..7dd0650180ac0cdc10a0cdceb847e0fce33fcda3 100644 (file)
@@ -845,14 +845,108 @@ namespace internal
         template <int spacedim>
         static
         void
-        renumber_mg_dofs (const std::vector<dealii::types::global_dof_index> &,
-                          const IndexSet &,
-                          DoFHandler<3,spacedim> &,
-                          const unsigned int               ,
-                          const bool                       )
+        renumber_mg_dofs (const std::vector<dealii::types::global_dof_index> &new_numbers,
+                          const IndexSet &indices,
+                          DoFHandler<3,spacedim> &dof_handler,
+                          const unsigned int level,
+                          const bool check_validity)
         {
-          // TODO
-          AssertThrow(false, ExcNotImplemented());
+          if (level>=dof_handler.get_tria().n_levels())
+            return;
+          for (typename std::vector<typename DoFHandler<3,spacedim>::MGVertexDoFs>::iterator i=dof_handler.mg_vertex_dofs.begin();
+               i!=dof_handler.mg_vertex_dofs.end(); ++i)
+            // if the present vertex lives on
+            // the present level
+            if ((i->get_coarsest_level() <= level) &&
+                (i->get_finest_level() >= level))
+              for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_vertex; ++d)
+                {
+                  dealii::types::global_dof_index idx =i->get_index (level, d);
+                  if (idx != DoFHandler<3>::invalid_dof_index)
+                    i->set_index (level, d,
+                                  ((indices.n_elements() == 0) ?
+                                   new_numbers[idx] :
+                                   new_numbers[indices.index_within_set(idx)]));
+
+                  if (check_validity)
+                    Assert(idx != DoFHandler<3>::invalid_dof_index, ExcInternalError ());
+                }
+
+          if (dof_handler.get_fe().dofs_per_line > 0 ||
+              dof_handler.get_fe().dofs_per_quad > 0)
+            {
+              // save user flags as they will be modified
+              std::vector<bool> user_flags;
+              dof_handler.get_tria().save_user_flags(user_flags);
+              const_cast<dealii::Triangulation<3,spacedim> &>(dof_handler.get_tria()).clear_user_flags ();
+
+              // flag all lines adjacent to cells of the current level, as
+              // those lines logically belong to the same level as the cell,
+              // at least for isotropic refinement
+              for (typename DoFHandler<3,spacedim>::level_cell_iterator cell = dof_handler.begin(level);
+                   cell != dof_handler.end(level); ++cell)
+                for (unsigned int line=0; line < GeometryInfo<3>::lines_per_cell; ++line)
+                  cell->line(line)->set_user_flag();
+
+              for (typename DoFHandler<3,spacedim>::cell_iterator cell = dof_handler.begin();
+                   cell != dof_handler.end(); ++cell)
+                for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
+                  if (cell->line(l)->user_flag_set())
+                    {
+                      for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_line; ++d)
+                        {
+                          dealii::types::global_dof_index idx = cell->line(l)->mg_dof_index(level, d);
+                          if (idx != DoFHandler<3>::invalid_dof_index)
+                            cell->line(l)->set_mg_dof_index (level, d, ((indices.n_elements() == 0) ?
+                                                                        new_numbers[idx] :
+                                                                        new_numbers[indices.index_within_set(idx)]));
+                          if (check_validity)
+                            Assert(idx != DoFHandler<3>::invalid_dof_index, ExcInternalError ());
+                        }
+                      cell->line(l)->clear_user_flag();
+                    }
+
+              // flag all quads adjacent to cells of the current level, as
+              // those quads logically belong to the same level as the cell,
+              // at least for isotropic refinement
+              for (typename DoFHandler<3,spacedim>::level_cell_iterator cell = dof_handler.begin(level);
+                   cell != dof_handler.end(level); ++cell)
+                for (unsigned int quad=0; quad < GeometryInfo<3>::quads_per_cell; ++quad)
+                  cell->quad(quad)->set_user_flag();
+
+              for (typename DoFHandler<3,spacedim>::cell_iterator cell = dof_handler.begin();
+                   cell != dof_handler.end(); ++cell)
+                for (unsigned int l=0; l<GeometryInfo<3>::quads_per_cell; ++l)
+                  if (cell->quad(l)->user_flag_set())
+                    {
+                      for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_quad; ++d)
+                        {
+                          dealii::types::global_dof_index idx = cell->quad(l)->mg_dof_index(level, d);
+                          if (idx != DoFHandler<1>::invalid_dof_index)
+                            cell->quad(l)->set_mg_dof_index (level, d, ((indices.n_elements() == 0) ?
+                                                                        new_numbers[idx] :
+                                                                        new_numbers[indices.index_within_set(idx)]));
+                          if (check_validity)
+                            Assert(idx != DoFHandler<3>::invalid_dof_index, ExcInternalError ());
+                        }
+                      cell->quad(l)->clear_user_flag();
+                    }
+
+              // finally, restore user flags
+              const_cast<dealii::Triangulation<3,spacedim> &>(dof_handler.get_tria()).load_user_flags (user_flags);
+            }
+
+          for (std::vector<types::global_dof_index>::iterator i=dof_handler.mg_levels[level]->dof_object.dofs.begin();
+               i!=dof_handler.mg_levels[level]->dof_object.dofs.end(); ++i)
+            {
+              if (*i != DoFHandler<3>::invalid_dof_index)
+                {
+                  Assert(*i<new_numbers.size(), ExcInternalError());
+                  *i = ((indices.n_elements() == 0) ?
+                        new_numbers[*i] :
+                        new_numbers[indices.index_within_set(*i)]);
+                }
+            }
         }
 
 

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.