]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix DoFTools::extract_constant_modes for FE_Q_DG0
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 4 Apr 2014 12:37:11 +0000 (12:37 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 4 Apr 2014 12:37:11 +0000 (12:37 +0000)
git-svn-id: https://svn.dealii.org/trunk@32714 0785d39b-7218-0410-832d-ea1e28bc413d

16 files changed:
deal.II/doc/news/changes.h
deal.II/include/deal.II/fe/fe.h
deal.II/include/deal.II/fe/fe_dgp.h
deal.II/include/deal.II/fe/fe_dgq.h
deal.II/include/deal.II/fe/fe_face.h
deal.II/include/deal.II/fe/fe_nedelec.h
deal.II/include/deal.II/fe/fe_q_base.h
deal.II/include/deal.II/fe/fe_q_dg0.h
deal.II/include/deal.II/fe/fe_q_hierarchical.h
deal.II/include/deal.II/fe/fe_raviart_thomas.h
deal.II/include/deal.II/fe/fe_system.h
deal.II/source/dofs/dof_tools.cc
tests/fe/element_constant_modes.cc
tests/fe/element_constant_modes.output
tests/mpi/extract_constant_modes_02.cc [new file with mode: 0644]
tests/mpi/extract_constant_modes_02.mpirun=4.output [new file with mode: 0644]

index c45afc16b7d3645ff4b145b0c821e35d1f17fbaf..a49dcacfe22503970b00c6b0ae0ff4b18d9e1c6a 100644 (file)
@@ -148,6 +148,14 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: DoFTools::extract_constant_modes now correctly identifies both
+  constant modes in the scalar element FE_Q_DG0, which has been realized by a
+  few modifications in how the constant modes propagate from the element to
+  the extract_constant_modes function.
+  <br>
+  (Martin Kronbichler, 2014/04/04)
+  </li>
+
   <li> Fixed: GridTools::laplace_transform had previously announced in
   the documentation that one can also set the location of interior points,
   but this was not in fact what was implemented. This has now been fixed:
index 84604ffc7e63f299dd5620d8f89538e0417de41f..1be926e630b3694c7133ce09a97ab4673f7f5b34 100644 (file)
@@ -1314,13 +1314,22 @@ public:
   block_mask (const ComponentMask &component_mask) const;
 
   /**
-   * Returns a list of constant modes of the element. The returns table has as
-   * many rows as there are components in the element and dofs_per_cell
-   * columns. To each component of the finite element, the row in the returned
-   * table contains a basis representation of the constant function 1 on the
-   * element.
-   */
-  virtual Table<2,bool> get_constant_modes () const;
+   * Returns a list of constant modes of the element. The number of rows in
+   * the resulting table depends on the elements in use. For standard
+   * elements, the table has as many rows as there are components in the
+   * element and dofs_per_cell columns. To each component of the finite
+   * element, the row in the returned table contains a basis representation of
+   * the constant function 1 on the element. However, there are some scalar
+   * elements where there is more than one constant mode, e.g. the element
+   * FE_Q_DG0.
+   *
+   * In order to match the constant modes to the actual components in the
+   * element, the returned data structure also returns a vector with as many
+   * components as there are constant modes on the element that contains the
+   * component number.
+   */
+  virtual std::pair<Table<2,bool>,std::vector<unsigned int> >
+  get_constant_modes () const;
 
   //@}
 
index 6343a7c637fdc102dc5c623fc09fdf4612df2a39..9fd6a5f04a7f87930a61970c484574343f65f95d 100644 (file)
@@ -331,7 +331,8 @@ public:
    * Returns a list of constant modes of the element. For this element, the
    * first entry is true, all other are false.
    */
-  virtual Table<2,bool> get_constant_modes () const;
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
 
 protected:
 
index 103c135c7e6c2257da58a277bcd5b85d21c21040..ec0c7f9e18abd3049e31bc6a43a83cd1d8ed60d3 100644 (file)
@@ -331,7 +331,8 @@ public:
    * Returns a list of constant modes of the element. For this element, it
    * simply returns one row with all entries set to true.
    */
-  virtual Table<2,bool> get_constant_modes () const;
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
 
   /**
    * Determine an estimate for the
index f6bd4606b5027dad2e3c254d8601fe20536da107..2ee366f50fb55bd77443d0eba8320d3848d3e538 100644 (file)
@@ -128,7 +128,8 @@ public:
    * Returns a list of constant modes of the element. For this element, it
    * simply returns one row with all entries set to true.
    */
-  virtual Table<2,bool> get_constant_modes () const;
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
 
 private:
   /**
@@ -230,7 +231,8 @@ public:
    * Returns a list of constant modes of the element. For this element, it
    * simply returns one row with all entries set to true.
    */
-  virtual Table<2,bool> get_constant_modes () const;
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
 
 protected:
   virtual
@@ -437,9 +439,10 @@ public:
    * Returns a list of constant modes of the element. For this element, the
    * first entry on each face is true, all other are false (as the constant
    * function is represented by the first base function of Legendre
-   * polynomials.
+   * polynomials).
    */
-  virtual Table<2,bool> get_constant_modes () const;
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
 
 private:
   /**
index 81dad103af2f725db6406de9bce67b5921cb2ba8..f09e5010bdb333cb1b71996b9b8a5ebd6c4c3980 100644 (file)
@@ -285,7 +285,8 @@ public:
   /**
    * Returns a list of constant modes of the element.
    */
-  virtual Table<2,bool> get_constant_modes () const;
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
 
   virtual std::size_t memory_consumption () const;
   virtual FiniteElement<dim> *clone() const;
index 6870f28fd27a8387dda487ae681cb8eb1d8a21ee..6bfb8a2ace2cabdd67e843670e866a970c8e02d4 100644 (file)
@@ -199,7 +199,8 @@ public:
    * Returns a list of constant modes of the element. For this element, the
    * list consists of true arguments for all components.
    */
-  virtual Table<2,bool> get_constant_modes () const;
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
 
   /**
    * @name Functions to support hp
index 6f0f339e01fc1d9a6fb7eecc666c831fdda63bf0..90b28f579acfb6cefe7e32fbc463ac1496490880 100644 (file)
@@ -312,6 +312,15 @@ public:
   virtual bool has_support_on_face (const unsigned int shape_index,
                                     const unsigned int face_index) const;
 
+  /**
+   * Returns a list of constant modes of the element. For this element, there
+   * are two constant modes despite the element is scalar: The first constant
+   * mode is all ones for the usual FE_Q basis and the second one only using
+   * the discontinuous part.
+   */
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
+
 protected:
   /**
    * @p clone function instead of a copy constructor.
index fd6602c5b11757ca88f15807f59227274e183e06..64d3b644cf97c1d6f1e6fc8618875c3d36321762 100644 (file)
@@ -386,7 +386,8 @@ public:
    * list consists of true arguments for the first vertex shape functions and
    * false for the remaining ones.
    */
-  virtual Table<2,bool> get_constant_modes () const;
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
 
 protected:
   /**
index bda93239da3c6ecc1662a9eca40d16d4ecd9acc6..01182aef014b9c49e754e6613b2a373474408ca1 100644 (file)
@@ -148,10 +148,11 @@ public:
     const VectorSlice<const std::vector<std::vector<double> > > &values) const;
 
   /**
-   * Returns a list of constant modes of the element. For this element, the
-   * list consists of true arguments for all components.
+   * Returns a list of constant modes of the element. This method is currently
+   * not correctly implemented because it returns ones for all components.
    */
-  virtual Table<2,bool> get_constant_modes () const;
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
 
   virtual std::size_t memory_consumption () const;
   virtual FiniteElement<dim> *clone() const;
index 7022cdd694e57a3940590bcdb7baeeece25506ec..04f97366889d742bde503983aa49640aec96e2a7 100644 (file)
@@ -476,7 +476,8 @@ public:
    * table contains a basis representation of the constant function 1 on the
    * element. Concatenates the constant modes of each base element.
    */
-  virtual Table<2,bool> get_constant_modes () const;
+  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
+  get_constant_modes () const;
 
   /**
    * @name Functions to support hp
index 3ea549933447ed05ab9bf7388230f450e8851844..1ca207eb0da10f288cfa5ac6005bdaefdde4bb09 100644 (file)
@@ -984,12 +984,6 @@ namespace DoFTools
     Assert (component_mask.represents_n_components(n_components),
             ExcDimensionMismatch(n_components,
                                  component_mask.size()));
-    std::vector<unsigned int> localized_component (n_components,
-                                                   numbers::invalid_unsigned_int);
-    unsigned int n_components_selected = 0;
-    for (unsigned int i=0; i<n_components; ++i)
-      if (component_mask[i] == true)
-        localized_component[i] = n_components_selected++;
 
     std::vector<unsigned char> dofs_by_component (dof_handler.n_locally_owned_dofs());
     internal::get_component_association (dof_handler, component_mask,
@@ -1008,17 +1002,37 @@ namespace DoFTools
       if (component_mask[dofs_by_component[i]])
         component_numbering[i] = count++;
 
-    // First count the number of dofs in the current component.
-    constant_modes.clear ();
-    constant_modes.resize (n_components_selected, std::vector<bool>(n_selected_dofs,
-                           false));
-
-    // Loop over all owned cells and ask the element for the constant modes
+    // get the element constant modes and find a translation table between
+    // index in the constant modes and the components.
+    //
+    // TODO: We might be able to extend this also for elements which do not
+    // have the same constant modes, but that is messy...
     const dealii::hp::FECollection<DH::dimension,DH::space_dimension>
       fe_collection (dof_handler.get_fe());
     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)
-      element_constant_modes.push_back(fe_collection[f].get_constant_modes());
+      {
+        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]].
+                push_back(std::make_pair(n_constant_modes++,i));
+        AssertDimension(element_constant_modes.back().n_rows(),
+                        element_constant_modes[0].n_rows());
+      }
+
+    // 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));
+
+    // Loop over all owned cells and ask the element for the constant modes
 
     typename DH::active_cell_iterator cell = dof_handler.begin_active(),
                                       endc = dof_handler.end();
@@ -1036,8 +1050,11 @@ namespace DoFTools
                   locally_owned_dofs.index_within_set(dof_indices[i]);
                 const unsigned int comp = dofs_by_component[loc_index];
                 if (component_mask[comp])
-                  constant_modes[localized_component[comp]][component_numbering[loc_index]] =
-                    element_constant_modes[cell->active_fe_index()](comp,i);
+                  for (unsigned int j=0; j<constant_mode_to_component_translation[comp].size(); ++j)
+                    constant_modes[constant_mode_to_component_translation[comp][j].first]
+                      [component_numbering[loc_index]] =
+                      element_constant_modes[cell->active_fe_index()]
+                      (constant_mode_to_component_translation[comp][j].second,i);
               }
         }
   }
index 4ae2eab1fa8898590da15a23c3a6ac4814711195..54e7d407c3ea2b502882455e7ecc665319ce80ff 100644 (file)
@@ -21,6 +21,7 @@
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_dgp.h>
 #include <deal.II/fe/fe_face.h>
+#include <deal.II/fe/fe_q_dg0.h>
 #include <deal.II/fe/fe_q_hierarchical.h>
 #include <deal.II/fe/fe_system.h>
 #include <fstream>
@@ -30,9 +31,9 @@
 template<int dim>
 void print_constant_modes(const FiniteElement<dim> &fe)
 {
-  Table<2,bool> constant_modes = fe.get_constant_modes();
   deallog << "Testing " << fe.get_name() << std::endl;
 
+  Table<2,bool> constant_modes = fe.get_constant_modes().first;
   for (unsigned int r=0; r<constant_modes.n_rows(); ++r)
     {
       for (unsigned int c=0; c<constant_modes.n_cols(); ++c)
@@ -56,6 +57,9 @@ void test()
   print_constant_modes(FE_FaceP<dim>(1));
   print_constant_modes(FESystem<dim>(FE_Q<dim>(1), 2, FE_Q<dim>(2), 1));
   print_constant_modes(FESystem<dim>(FE_DGP<dim>(1), 1, FE_Q_iso_Q1<dim>(2), 1));
+  print_constant_modes(FE_Q_DG0<dim>(1));
+  print_constant_modes(FESystem<dim>(FE_Q_DG0<dim>(2), 1, FE_Q<dim>(1), 2));
+  print_constant_modes(FESystem<dim>(FE_Q<dim>(1), 2, FE_Q_DG0<dim>(1), 2));
 }
 
 template <>
index 428dbf5c6e48fbbea700df81afc2b9048c0d922c..076019878b56462d6adc5d5605d22b5e56877cfe 100644 (file)
@@ -44,6 +44,24 @@ DEAL::Testing FESystem<2>[FE_DGP<2>(1)-FE_Q_iso_Q1<2>(2)]
 DEAL::0 0 0 0 0 0 0 0 1 0 0 0 
 DEAL::1 1 1 1 1 1 1 1 0 0 0 1 
 DEAL::
+DEAL::Testing FE_Q_DG0<2>(1)
+DEAL::1 1 1 1 0 
+DEAL::0 0 0 0 1 
+DEAL::
+DEAL::Testing FESystem<2>[FE_Q_DG0<2>(2)-FE_Q<2>(1)^2]
+DEAL::1 0 0 1 0 0 1 0 0 1 0 0 1 1 1 1 1 0 
+DEAL::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 
+DEAL::0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 
+DEAL::0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 
+DEAL::
+DEAL::Testing FESystem<2>[FE_Q<2>(1)^2-FE_Q_DG0<2>(1)^2]
+DEAL::1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 
+DEAL::0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 
+DEAL::0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 
+DEAL::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 
+DEAL::0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 
+DEAL::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 
+DEAL::
 DEAL::Testing FE_Q<3>(1)
 DEAL::1 1 1 1 1 1 1 1 
 DEAL::
@@ -77,3 +95,21 @@ DEAL::Testing FESystem<3>[FE_DGP<3>(1)-FE_Q_iso_Q1<3>(2)]
 DEAL::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 
 DEAL::1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 
 DEAL::
+DEAL::Testing FE_Q_DG0<3>(1)
+DEAL::1 1 1 1 1 1 1 1 0 
+DEAL::0 0 0 0 0 0 0 0 1 
+DEAL::
+DEAL::Testing FESystem<3>[FE_Q_DG0<3>(2)-FE_Q<3>(1)^2]
+DEAL::1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 
+DEAL::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 
+DEAL::0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::
+DEAL::Testing FESystem<3>[FE_Q<3>(1)^2-FE_Q_DG0<3>(1)^2]
+DEAL::1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 
+DEAL::0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 
+DEAL::0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 
+DEAL::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 
+DEAL::0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 
+DEAL::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 
+DEAL::
diff --git a/tests/mpi/extract_constant_modes_02.cc b/tests/mpi/extract_constant_modes_02.cc
new file mode 100644 (file)
index 0000000..d03cd60
--- /dev/null
@@ -0,0 +1,104 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2009 - 2013 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test DoFTools::extract_constant_modes with FE_Q_DG0
+
+#include "../tests.h"
+#include "coarse_grid_common.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <deal.II/fe/fe_q_dg0.h>
+
+#include <fstream>
+
+
+template<int dim>
+void test()
+{
+  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] = 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));
+
+  FE_Q_DG0<dim> fe(1);
+  DoFHandler<dim> dofh(tr);
+  dofh.distribute_dofs (fe);
+
+  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(fe.n_components(), 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)
+    {
+      std::ofstream logfile("output");
+      deallog.attach(logfile);
+      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+
+      deallog.push("2d");
+      test<2>();
+      deallog.pop();
+    }
+  else
+    {
+      deallog.push("2d");
+      test<2>();
+      deallog.pop();
+    }
+}
diff --git a/tests/mpi/extract_constant_modes_02.mpirun=4.output b/tests/mpi/extract_constant_modes_02.mpirun=4.output
new file mode 100644 (file)
index 0000000..63cce8e
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL:0:2d::Total dofs=14
+DEAL:0:2d::1 1 1 1 0 
+DEAL:0:2d::0 0 0 0 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.