]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement DoFTools::count_dofs_per_component also for hp::DoFhandler.
authorgoll <goll@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Jul 2011 16:14:35 +0000 (16:14 +0000)
committergoll <goll@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Jul 2011 16:14:35 +0000 (16:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@23924 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/dofs/dof_tools.h
deal.II/source/dofs/dof_tools.cc
deal.II/source/dofs/dof_tools.inst.in
tests/bits/count_dofs_per_component_hp.cc [new file with mode: 0644]
tests/bits/count_dofs_per_component_hp/cmp/generic [new file with mode: 0644]

index 8a0a0fe2ba751a54029a9bfd3507b2fe155e286a..1dbab530f86f31c1c6d0bf37fecf6f6ed427e78d 100644 (file)
@@ -216,6 +216,11 @@ should be fixed now.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> New: The function DoFTools::count_dofs_per_component now also works
+for objects of type hp::DoFHandler.
+<br>
+(Christian Goll, Wolfgang Bangerth 2011/07/06)
+
 <li> Fixed: Under some circumstances, the Threads::Thread::join() could only be
 called once and would generate a system exception when called a second time.
 Since it is often useful to not track whether this function had already been
index 26585a21930f1747d73ae24574a75bdca67542e4..2f9b90f2b3a496a771a09e388b69a11b7ccc69b7 100644 (file)
@@ -1580,9 +1580,9 @@ namespace DoFTools
                                    * targetted by target_components
                                    * are left untouched.
                                    */
-  template <int dim, int spacedim>
+  template <class DH>
   void
-  count_dofs_per_component (const DoFHandler<dim,spacedim>&     dof_handler,
+  count_dofs_per_component (const D&     dof_handler,
                            std::vector<unsigned int>& dofs_per_component,
                            const bool vector_valued_once = false,
                            std::vector<unsigned int>  target_component
@@ -2100,7 +2100,7 @@ namespace DoFTools
   std::vector<Table<2,Coupling> >
   dof_couplings_from_component_couplings (const hp::FECollection<dim,spacedim> &fe,
                                          const Table<2,Coupling> &component_couplings);
-    
+
                                   /**
                                    * Exception
                                    */
index 56003aa3f4f3d3a8ad4060368f3f5b069223620e..335931698db4a8a45a340298de9a9a45a920ea8f 100644 (file)
@@ -51,7 +51,7 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace DoFTools
 {
-  
+
   template <class DH, class SparsityPattern>
   void
   make_sparsity_pattern (const DH               &dof,
@@ -732,14 +732,14 @@ namespace DoFTools
     return return_value;
   }
 
-  
+
 
 
   namespace internal
   {
-    namespace 
+    namespace
     {
-      
+
                                       // implementation of the same function in
                                       // namespace DoFTools for non-hp
                                       // DoFHandlers
@@ -972,7 +972,7 @@ namespace DoFTools
                }
             }
       }
-        
+
 
                                       // implementation of the same function in
                                       // namespace DoFTools for non-hp
@@ -1166,14 +1166,14 @@ namespace DoFTools
          }
       }
     }
-    
+
   }
 
 
 
 
   template <class DH, class SparsityPattern>
-  void  
+  void
   make_flux_sparsity_pattern (const DH                &dof,
                              SparsityPattern         &sparsity,
                              const Table<2,Coupling> &int_mask,
@@ -1199,7 +1199,7 @@ namespace DoFTools
     Assert (flux_mask.n_cols() == n_comp,
            ExcDimensionMismatch (flux_mask.n_cols(), n_comp));
 
-  
+
                                     // Clear user flags because we will
                                     // need them. But first we save
                                     // them and make sure that we
@@ -1214,7 +1214,7 @@ namespace DoFTools
 
     internal::make_flux_sparsity_pattern (dof, sparsity,
                                          int_mask, flux_mask);
-  
+
                                     // finally restore the user flags
     const_cast<Triangulation<DH::dimension,DH::space_dimension> &>(dof.get_tria()).load_user_flags(user_flags);
   }
@@ -2710,7 +2710,7 @@ namespace DoFTools
                              ==
                              FiniteElementDomination::no_requirements)
                            continue;
-                         
+
                                                           // Same procedure as for the
                                                           // mother cell. Extract the face
                                                           // DoFs from the cell DoFs.
@@ -3245,7 +3245,7 @@ namespace DoFTools
                                                           // nothing to do here
                          break;
                        }
-                       
+
                        default:
                                                               // we shouldn't get
                                                               // here
@@ -4085,7 +4085,7 @@ namespace DoFTools
 
   template <int dim, int spacedim>
   void
-  
+
   extract_hanging_node_dofs (const DoFHandler<dim,spacedim> &dof_handler,
                             std::vector<bool>              &selected_dofs)
   {
@@ -4470,6 +4470,10 @@ namespace DoFTools
 
   namespace internal
   {
+//TODO: why is this function so complicated? It would be nice to have
+// comments that explain why we can't just loop over all components
+// and count the entries in dofs_by_component that have this
+// component's index
     template <int dim, int spacedim>
     void
     resolve_components (const FiniteElement<dim,spacedim>&fe,
@@ -4509,18 +4513,78 @@ namespace DoFTools
            }
        }
     }
+
+    template <int dim, int spacedim>
+    void
+    resolve_components (const hp::FECollection<dim,spacedim>&fe_collection,
+                        const std::vector<unsigned char> &dofs_by_component,
+                        const std::vector<unsigned int>  &target_component,
+                        const bool                        only_once,
+                        std::vector<unsigned int>        &dofs_per_component,
+                        unsigned int                     &component)
+    {
+      // assert that all elements in the collection have the same
+      // structure (base elements and multiplicity, components per base
+      // element) and then simply call the function above
+      for (unsigned int fe=1; fe<fe_collection.size(); ++fe)
+      {
+        Assert (fe_collection[fe].n_components() == fe_collection[0].n_components(),
+                ExcNotImplemented());
+        Assert (fe_collection[fe].n_base_elements() == fe_collection[0].n_base_elements(),
+                ExcNotImplemented());
+        for (unsigned int b=0;b<fe_collection[0].n_base_elements();++b)
+        {
+          Assert (fe_collection[fe].base_element(b).n_components() == fe_collection[0].base_element(b).n_components(),
+                  ExcNotImplemented());
+          Assert (fe_collection[fe].base_element(b).n_base_elements() == fe_collection[0].base_element(b).n_base_elements(),
+                  ExcNotImplemented());
+        }
+      }
+
+      resolve_components (fe_collection[0], dofs_by_component,
+                          target_component, only_once, dofs_per_component,
+                          component);
+    }
   }
 
 
-  template <int dim, int spacedim>
+  namespace internal
+  {
+    namespace {
+      /**
+       * Return true if the given element is primitive.
+       */
+      template <int dim, int spacedim>
+      bool all_elements_are_primitive (const FiniteElement<dim,spacedim> &fe)
+      {
+        return fe.is_primitive();
+      }
+
+
+      /**
+       * Return true if each element of the given element collection is primitive.
+       */
+      template <int dim, int spacedim>
+      bool all_elements_are_primitive (const dealii::hp::FECollection<dim,spacedim> &fe_collection)
+      {
+        for (unsigned int i=0; i<fe_collection.size(); ++i)
+          if (fe_collection[i].is_primitive() == false)
+            return false;
+
+        return true;
+      }
+    }
+  }
+
+  template <class DH>
   void
   count_dofs_per_component (
-    const DoFHandler<dim,spacedim>&     dof_handler,
+    const D&     dof_handler,
     std::vector<unsigned int>& dofs_per_component,
     bool only_once,
     std::vector<unsigned int>  target_component)
   {
-    const FiniteElement<dim,spacedim>& fe = dof_handler.get_fe();
+    const unsigned int n_components = dof_handler.get_fe().n_components();
 
     std::fill (dofs_per_component.begin(), dofs_per_component.end(), 0U);
 
@@ -4529,21 +4593,20 @@ namespace DoFTools
                                     // vector as identity.
     if (target_component.size()==0)
       {
-       target_component.resize(fe.n_components());
-       for (unsigned int i=0; i<fe.n_components(); ++i)
+       target_component.resize(n_components);
+       for (unsigned int i=0; i<n_components; ++i)
          target_component[i] = i;
       }
     else
-      Assert (target_component.size()==fe.n_components(),
+      Assert (target_component.size()==n_components,
              ExcDimensionMismatch(target_component.size(),
-                                  fe.n_components()));
+                                  n_components));
 
 
     const unsigned int max_component
       = *std::max_element (target_component.begin(),
                           target_component.end());
     const unsigned int n_target_components = max_component + 1;
-    const unsigned int n_components = fe.n_components();
 
     AssertDimension (dofs_per_component.size(), n_target_components);
 
@@ -4567,7 +4630,8 @@ namespace DoFTools
 
                                     // next count what we got
     unsigned int component = 0;
-    internal::resolve_components(fe, dofs_by_component, target_component,
+    internal::resolve_components(dof_handler.get_fe(),
+                                 dofs_by_component, target_component,
                                 only_once, dofs_per_component, component);
     Assert (n_components == component, ExcInternalError());
 
@@ -4575,7 +4639,7 @@ namespace DoFTools
                                     // only valid if the finite element
                                     // is actually primitive, so
                                     // exclude other elements from this
-    Assert (!dof_handler.get_fe().is_primitive()
+    Assert ((internal::all_elements_are_primitive(dof_handler.get_fe()) == false)
            ||
            (std::accumulate (dofs_per_component.begin(),
                              dofs_per_component.end(), 0U)
@@ -4585,8 +4649,8 @@ namespace DoFTools
                                     // reduce information from all CPUs
 #ifdef DEAL_II_USE_P4EST
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-    if (const parallel::distributed::Triangulation<dim> * tria
-       = (dynamic_cast<const parallel::distributed::Triangulation<dim>*>
+    if (const parallel::distributed::Triangulation<dim,spacedim> * tria
+       = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
           (&dof_handler.get_tria())))
       {
        std::vector<unsigned int> local_dof_count = dofs_per_component;
@@ -4601,7 +4665,7 @@ namespace DoFTools
 
 
   template <int dim, int spacedim>
-  void  
+  void
   count_dofs_per_block (const DoFHandler<dim,spacedim>& dof_handler,
                        std::vector<unsigned int> &dofs_per_block,
                        std::vector<unsigned int>  target_block)
index ddabf8f8c28a50d35cfbc2311e80732db0cfbfd6..28287e1cf924f88742c2ab84003b38a6bd4a993b 100644 (file)
@@ -487,20 +487,36 @@ DoFTools::count_dofs_with_subdomain_association<MGDoFHandler<deal_II_dimension>
 
 template
 void
-DoFTools::count_dofs_per_component<deal_II_dimension> (
+DoFTools::count_dofs_per_component<DoFHandler<deal_II_dimension> > (
   const DoFHandler<deal_II_dimension>&,
   std::vector<unsigned int>&, bool, std::vector<unsigned int>);
 
-#if deal_II_dimension < 3
-template void DoFTools::extract_level_dofs<deal_II_dimension, deal_II_dimension+1>
-(const unsigned int level, const MGDoFHandler<deal_II_dimension, deal_II_dimension+1>&,
- const std::vector<bool>&, std::vector<bool>&, bool);
+template
+void
+DoFTools::count_dofs_per_component<hp::DoFHandler<deal_II_dimension> > (
+  const hp::DoFHandler<deal_II_dimension>&,
+  std::vector<unsigned int>&, bool, std::vector<unsigned int>);
+
 
+#if deal_II_dimension < 3
 template
 void
-DoFTools::count_dofs_per_component<deal_II_dimension, deal_II_dimension+1> (
+DoFTools::count_dofs_per_component<DoFHandler<deal_II_dimension, deal_II_dimension+1> > (
   const DoFHandler<deal_II_dimension, deal_II_dimension+1>&,
   std::vector<unsigned int>&, bool, std::vector<unsigned int>);
+
+template
+void
+DoFTools::count_dofs_per_component<hp::DoFHandler<deal_II_dimension, deal_II_dimension+1> > (
+  const hp::DoFHandler<deal_II_dimension, deal_II_dimension+1>&,
+  std::vector<unsigned int>&, bool, std::vector<unsigned int>);
+  
+template 
+void 
+DoFTools::extract_level_dofs<deal_II_dimension, deal_II_dimension+1>
+(const unsigned int level, const MGDoFHandler<deal_II_dimension, deal_II_dimension+1>&,
+ const std::vector<bool>&, std::vector<bool>&, bool);
+  
 #endif
 
 template
diff --git a/tests/bits/count_dofs_per_component_hp.cc b/tests/bits/count_dofs_per_component_hp.cc
new file mode 100644 (file)
index 0000000..184f139
--- /dev/null
@@ -0,0 +1,83 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2003, 2004, 2005, 2009 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------------
+
+#include "../tests.h"
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <hp/fe_collection.h>
+#include <hp/dof_handler.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/fe_system.h>
+#include <dofs/dof_tools.h>
+
+// check
+//   DoFTools::
+//   count_dofs_per_component (...);
+// for the hp case
+
+
+using namespace std;
+
+template <int dim>
+void test ()
+{
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation);
+  triangulation.refine_global(3);
+
+  //define DoFhandler and FEs
+  FE_Q<dim> u(2);
+  FE_Q<dim> p(1);
+
+  FESystem<dim> fe_system(u, 2, p, 1);
+
+  hp::FECollection<dim> fe_collection;
+  fe_collection.push_back(fe_system);
+
+  hp::DoFHandler<dim> hp_dof_handler(triangulation);
+  DoFHandler<dim> dof_handler(triangulation);
+
+  //distribute dofs
+  hp_dof_handler.distribute_dofs(fe_collection);
+  dof_handler.distribute_dofs(fe_system);
+
+  //count dofs per component and show them on the screen
+  std::vector<unsigned int> dofs_per_component(3,0);
+  std::vector<unsigned int> dofs_per_component_hp(3,0);
+  DoFTools::count_dofs_per_component(dof_handler, dofs_per_component);
+  DoFTools::count_dofs_per_component(hp_dof_handler, dofs_per_component_hp);
+
+  for (unsigned int i=0; i<3; i++)
+    {
+      deallog <<"DoFs in the " <<i<<". component for classical FE: "<< dofs_per_component.at(i)<< std::endl;
+      deallog <<"DoFs in the " <<i<<". component for hp FE: "<< dofs_per_component_hp.at(i)<< std::endl;
+
+      Assert (dofs_per_component.at(i) == dofs_per_component_hp.at(i),
+             ExcInternalError());
+    }
+
+}
+
+int main () 
+{
+  std::ofstream logfile("count_dofs_per_component_hp/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<1> ();
+  test<2> ();
+  test<3> ();
+  return 0;
+}
diff --git a/tests/bits/count_dofs_per_component_hp/cmp/generic b/tests/bits/count_dofs_per_component_hp/cmp/generic
new file mode 100644 (file)
index 0000000..d51c947
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL::DoFs in the 0. component for classical FE: 17
+DEAL::DoFs in the 0. component for hp FE: 17
+DEAL::DoFs in the 1. component for classical FE: 17
+DEAL::DoFs in the 1. component for hp FE: 17
+DEAL::DoFs in the 2. component for classical FE: 9
+DEAL::DoFs in the 2. component for hp FE: 9
+DEAL::DoFs in the 0. component for classical FE: 289
+DEAL::DoFs in the 0. component for hp FE: 289
+DEAL::DoFs in the 1. component for classical FE: 289
+DEAL::DoFs in the 1. component for hp FE: 289
+DEAL::DoFs in the 2. component for classical FE: 81
+DEAL::DoFs in the 2. component for hp FE: 81
+DEAL::DoFs in the 0. component for classical FE: 4913
+DEAL::DoFs in the 0. component for hp FE: 4913
+DEAL::DoFs in the 1. component for classical FE: 4913
+DEAL::DoFs in the 1. component for hp FE: 4913
+DEAL::DoFs in the 2. component for classical FE: 729
+DEAL::DoFs in the 2. component for hp FE: 729

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.