]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Enable FE_Nothing in DoFTools::extract_constant_modes 13313/head
authorSebastian Proell <sebastian.proell@tum.de>
Mon, 31 Jan 2022 17:13:58 +0000 (18:13 +0100)
committerSebastian Proell <sebastian.proell@tum.de>
Mon, 7 Feb 2022 08:27:43 +0000 (09:27 +0100)
doc/news/changes/minor/20220207Proell [new file with mode: 0644]
include/deal.II/fe/fe_nothing.h
source/dofs/dof_tools.cc
source/fe/fe_nothing.cc
tests/mpi/extract_constant_modes_03.cc [new file with mode: 0644]
tests/mpi/extract_constant_modes_03.mpirun=2.with_p4est=true.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20220207Proell b/doc/news/changes/minor/20220207Proell
new file mode 100644 (file)
index 0000000..0b3a1de
--- /dev/null
@@ -0,0 +1,5 @@
+New: The FE_Nothing class can now return its (empty) constant modes.
+DoFTools::extract_constant_modes can now work with DoFHandlers that
+contain FE_Nothing.
+<br>
+(Sebastian Proell, 2022/02/07)
index 2df7f0b5a52a4108df6cadb7b2d4e5d36d7de69d..ec229d6dc1d4dbc9858f6f42b69c326d2846fae0 100644 (file)
@@ -326,6 +326,15 @@ public:
     FullMatrix<double> &                interpolation_matrix,
     const unsigned int                  face_no = 0) const override;
 
+  /**
+   * Return a list of constant modes of the element.
+   *
+   * Since the current finite element has no degrees of freedom, the returned
+   * list is necessarily empty.
+   */
+  virtual std::pair<Table<2, bool>, std::vector<unsigned int>>
+  get_constant_modes() const override;
+
   /**
    * @return true if the FE dominates any other.
    */
index 1acf4d5602d531d7a50a68f5b1237a5e42dff9bd..ba7379924e5e1df2e0d55519d8cb1f381d1ff024 100644 (file)
@@ -1249,26 +1249,47 @@ namespace DoFTools
       dof_handler.get_fe_collection();
     std::vector<Table<2, bool>> element_constant_modes;
     std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
-                 constant_mode_to_component_translation(n_components);
-    unsigned int n_constant_modes = 0;
-    for (unsigned int f = 0; f < fe_collection.size(); ++f)
-      {
-        std::pair<Table<2, bool>, std::vector<unsigned int>> data =
-          fe_collection[f].get_constant_modes();
-        element_constant_modes.push_back(data.first);
-        if (f == 0)
-          for (unsigned int i = 0; i < data.second.size(); ++i)
-            if (component_mask[data.second[i]])
-              constant_mode_to_component_translation[data.second[i]]
-                .emplace_back(n_constant_modes++, i);
-        AssertDimension(element_constant_modes.back().n_rows(),
-                        element_constant_modes[0].n_rows());
-      }
+      constant_mode_to_component_translation(n_components);
+    {
+      unsigned int n_constant_modes              = 0;
+      int          first_non_empty_constant_mode = -1;
+      for (unsigned int f = 0; f < fe_collection.size(); ++f)
+        {
+          std::pair<Table<2, bool>, std::vector<unsigned int>> data =
+            fe_collection[f].get_constant_modes();
 
-    // First count the number of dofs in the current component.
-    constant_modes.clear();
-    constant_modes.resize(n_constant_modes,
-                          std::vector<bool>(n_selected_dofs, false));
+          // Store the index of the current element if it is the first that has
+          // non-empty constant modes.
+          if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0)
+            {
+              first_non_empty_constant_mode = f;
+              // This is the first non-empty constant mode, so we figure out the
+              // translation between index in the constant modes and the
+              // components
+              for (unsigned int i = 0; i < data.second.size(); ++i)
+                if (component_mask[data.second[i]])
+                  constant_mode_to_component_translation[data.second[i]]
+                    .emplace_back(n_constant_modes++, i);
+            }
+
+          // Add the constant modes of this element to the list and assert that
+          // there are as many constant modes as for the other elements (or zero
+          // constant modes).
+          element_constant_modes.push_back(data.first);
+          Assert(
+            element_constant_modes.back().n_rows() == 0 ||
+              element_constant_modes.back().n_rows() ==
+                element_constant_modes[first_non_empty_constant_mode].n_rows(),
+            ExcInternalError());
+        }
+      AssertIndexRange(first_non_empty_constant_mode, fe_collection.size());
+
+      // Now we know the number of constant modes and resize the return vector
+      // accordingly
+      constant_modes.clear();
+      constant_modes.resize(n_constant_modes,
+                            std::vector<bool>(n_selected_dofs, false));
+    }
 
     // Loop over all owned cells and ask the element for the constant modes
     std::vector<types::global_dof_index> dof_indices;
index 40db67be56547b517b22407645a3c0b6dfaec0a1..a1742e80b4184fcb690aa11d1eebdd07db74d0e1 100644 (file)
@@ -314,6 +314,7 @@ FE_Nothing<dim, spacedim>::get_face_interpolation_matrix(
 }
 
 
+
 template <int dim, int spacedim>
 void
 FE_Nothing<dim, spacedim>::get_subface_interpolation_matrix(
@@ -334,6 +335,16 @@ FE_Nothing<dim, spacedim>::get_subface_interpolation_matrix(
 
 
 
+template <int dim, int spacedim>
+std::pair<Table<2, bool>, std::vector<unsigned int>>
+FE_Nothing<dim, spacedim>::get_constant_modes() const
+{
+  // since this element has no dofs, there are no constant modes
+  return {Table<2, bool>{}, std::vector<unsigned int>{}};
+}
+
+
+
 // explicit instantiations
 #include "fe_nothing.inst"
 
diff --git a/tests/mpi/extract_constant_modes_03.cc b/tests/mpi/extract_constant_modes_03.cc
new file mode 100644 (file)
index 0000000..144eab2
--- /dev/null
@@ -0,0 +1,129 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test DoFTools::extract_constant_modes with FE_Nothing
+
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test(unsigned int fe_nothing_index)
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  parallel::distributed::Triangulation<dim> tr(MPI_COMM_WORLD);
+
+  std::vector<unsigned int> sub(2);
+  sub[0] = 2 * Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  sub[1] = 1;
+  GridGenerator::subdivided_hyper_rectangle(
+    static_cast<Triangulation<dim> &>(tr), sub, Point<2>(0, 0), Point<2>(1, 1));
+
+  DoFHandler<dim> dofh(tr);
+
+  {
+    // set cells to use FE_Q and FE_Nothing alternately
+    unsigned int last_index = 1;
+    for (const auto &cell :
+         dofh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+      cell->set_active_fe_index(last_index = 1 - last_index);
+  }
+
+  if (fe_nothing_index == 1)
+    dofh.distribute_dofs(
+      hp::FECollection<dim>(FE_Q<dim>(1), FE_Nothing<dim>(1)));
+  else
+    dofh.distribute_dofs(
+      hp::FECollection<dim>(FE_Nothing<dim>(1), FE_Q<dim>(1)));
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    deallog << "Total dofs=" << dofh.n_dofs() << std::endl;
+
+  // extract constant modes and print
+
+  if (myid == 0)
+    {
+      std::vector<bool> mask(1, true);
+
+      std::vector<std::vector<bool>> constant_modes;
+      DoFTools::extract_constant_modes(dofh, mask, constant_modes);
+
+      for (unsigned int i = 0; i < constant_modes.size(); ++i)
+        {
+          for (unsigned int j = 0; j < constant_modes[i].size(); ++j)
+            deallog << (constant_modes[i][j] ? '1' : '0') << ' ';
+          deallog << std::endl;
+        }
+    }
+}
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  deallog.push(Utilities::int_to_string(myid));
+
+
+  if (myid == 0)
+    {
+      initlog();
+
+      deallog.push("FE_Nothing second");
+      test<2>(1);
+      deallog.pop();
+
+      deallog.push("FE_Nothing first");
+      test<2>(0);
+      deallog.pop();
+    }
+  else
+    {
+      deallog.push("FE_Nothing second");
+      test<2>(1);
+      deallog.pop();
+
+      deallog.push("FE_Nothing first");
+      test<2>(0);
+      deallog.pop();
+    }
+}
diff --git a/tests/mpi/extract_constant_modes_03.mpirun=2.with_p4est=true.output b/tests/mpi/extract_constant_modes_03.mpirun=2.with_p4est=true.output
new file mode 100644 (file)
index 0000000..6b81914
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:0:FE_Nothing second::Total dofs=8
+DEAL:0:FE_Nothing second::1 1 1 1 
+DEAL:0:FE_Nothing first::Total dofs=8
+DEAL:0:FE_Nothing first::1 1 1 1 

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.