]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
multigrid tests improved
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 17 Mar 2010 17:44:27 +0000 (17:44 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 17 Mar 2010 17:44:27 +0000 (17:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@20841 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/multigrid/mg_constraints.h [new file with mode: 0644]
deal.II/deal.II/include/multigrid/mg_dof_accessor.h
deal.II/deal.II/source/multigrid/mg_tools.cc
deal.II/examples/step-39/step-39.cc

diff --git a/deal.II/deal.II/include/multigrid/mg_constraints.h b/deal.II/deal.II/include/multigrid/mg_constraints.h
new file mode 100644 (file)
index 0000000..dd009b9
--- /dev/null
@@ -0,0 +1,111 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2010 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__mg_constraints_h
+#define __deal2__mg_constraints_h
+
+#include <base/config.h>
+#include <base/subscriptor.h>
+
+#include <vector>
+#include <set>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * Collection of boundary constraints and refinement edge constraints
+ * for level vectors.
+ */
+class MGConstraints : public Subscriptor
+{
+  public:
+                                    /**
+                                     * Fill the internal data
+                                     * structures with values
+                                     * extracted from the dof
+                                     * handler.
+                                     *
+                                     * This function leaves
+                                     * #boundary_indices empty, since
+                                     * no boundary values are
+                                     * provided.
+                                     */
+    template <int dim, int spacedim>
+    void initialize(const MGDoFHandler<dim,spacedim>& dof);
+    
+                                    /**
+                                     * Fill the internal data
+                                     * structures with values
+                                     * extracted from the dof
+                                     * handler, applying the boundary
+                                     * values provided.
+                                     */
+    template <int dim, int spacedim>
+    void initialize(const MGDoFHandler<dim,spacedim>& dof,
+                   const typename FunctionMap<dim>::type& function_map,
+                   const std::vector<bool>& component_mask = std::vector<bool>());
+
+                                    /**
+                                     * Reset the data structures.
+                                     */
+    void clear();
+
+                                    /**
+                                     * Determine whether a dof index
+                                     * is subject to a boundary
+                                     * constraint.
+                                     */
+    bool at_boundary(unsigned int level, unsigned int index) const;
+    
+                                    /**
+                                     * Determine whether a dof index
+                                     * is at the refinement edge.
+                                     */
+    bool at_refinement_edge(unsigned int level, unsigned int index) const;
+    
+  private:
+                                    /**
+                                     * The indices of boundary dofs
+                                     * for each level.
+                                     */
+    std::vector<std::set<unsigned int> > boundary_indices;
+                                    /**
+                                     * The degrees of freedom on the
+                                     * refinement edge between a
+                                     * level and coarser cells.
+                                     */
+    std::vector<std::vector<bool> > refinement_edge_indices;
+    
+                                    /**
+                                     * The degrees of freedom on the
+                                     * refinement edge between a
+                                     * level and coarser cells, which
+                                     * are also on the boundary.
+                                     *
+                                     * This is a subset of
+                                     * #refinement_edge_indices.
+                                     */
+    std::vector<std::vector<bool> > refinement_edge_boundary_indices;
+};
+
+
+bool
+MGConstraints::at_boundary(unsigned int level, unsigned int index)
+{
+  AssertIndexRange(level, boundary_indices.size());
+  AssertIndexRange(level, boundary_indices.size());  
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 948ffc2d284320f9c7a9961278d81dc79ef5d75b..01a42d0fa0511526157695d9a34474f2db6d9af2 100644 (file)
@@ -289,8 +289,11 @@ class MGDoFAccessor : public internal::MGDoFAccessor::Inheritance<structdim,dim,
     unsigned int mg_vertex_dof_index (const int level,
                                      const unsigned int vertex,
                                      const unsigned int i) const;
-
+//TODO: This function is broken for faces
                                     /**
+                                     * @todo This function is broken
+                                     * for faces.
+                                     *
                                      * Return the index of the @p ith
                                      * degree of freedom of this line
                                      * on the level this line lives
index 805bc755978565f76e32b7dcd5cadd465aa949c8..a5673cdb2e9d30f2468e3136d76395071560839f 100644 (file)
@@ -1127,148 +1127,187 @@ MGTools::make_boundary_list(
   if (function_map.size() == 0)
     return;
 
-  const unsigned int n_components = DoFTools::n_components(dof);
   const unsigned int n_levels = dof.get_tria().n_levels();
+
+  
+  
+  const unsigned int n_components = DoFTools::n_components(dof);
   const bool          fe_is_system = (n_components != 1);
 
   AssertDimension (boundary_indices.size(), n_levels);
 
+  std::vector<unsigned int> local_dofs;
+  
+                                  // First, deal with the simpler
+                                  // case when we have to identify
+                                  // all boundary dofs
+  if (component_mask_.size() == 0)
+    {
+      typename MGDoFHandler<dim,spacedim>::cell_iterator
+       cell = dof.begin(),
+       endc = dof.end();
+      for (; cell!=endc; ++cell)
+       {
+         const FiniteElement<dim> &fe = cell->get_fe();
+         const unsigned int level = cell->level();
+         local_dofs.resize(fe.dofs_per_face);
+         
+         for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+              ++face_no)
+           {
+             const typename MGDoFHandler<dim,spacedim>::face_iterator
+               face = cell->face(face_no);
+             face->get_mg_dof_indices(level, local_dofs);
+             const unsigned char bi = face->boundary_indicator();
+                                              // Face is listed in
+                                              // boundary map
+             if (function_map.find(bi) != function_map.end())
+               {
+                 for (unsigned int i=0;i<fe.dofs_per_face;++i)
+                   boundary_indices[level].insert(i);
+               }
+           }
+       }
+    }
+  else
+    {
+      
 //    for (typename FunctionMap<dim>::type::const_iterator i=function_map.begin();
 //         i!=function_map.end(); ++i)
 //      Assert (n_components == i->second->n_components,
 //           ExcInvalidFE());
-
-                                   // set the component mask to either
-                                   // the original value or a vector
-                                   // of @p{true}s
-  const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
-                                          std::vector<bool> (n_components, true) :
-                                          component_mask_);
+      
+                                      // set the component mask to either
+                                      // the original value or a vector
+                                      // of @p{true}s
+      const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
+                                             std::vector<bool> (n_components, true) :
+                                             component_mask_);
 //   Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
 //        ExcComponentMismatch());
-
-                                   // field to store the indices
-  std::vector<unsigned int> face_dofs;
-  face_dofs.reserve (DoFTools::max_dofs_per_face(dof));
-  std::fill (face_dofs.begin (), face_dofs.end (),
-            DoFHandler<dim,spacedim>::invalid_dof_index);
-
-  typename MGDoFHandler<dim,spacedim>::cell_iterator
-    cell = dof.begin(),
-    endc = dof.end();
-  for (; cell!=endc; ++cell)
-    for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
-         ++face_no)
-      {
-        const FiniteElement<dim> &fe = cell->get_fe();
-        const unsigned int level = cell->level();
-
-                                         // we can presently deal only with
-                                         // primitive elements for boundary
-                                         // values. this does not preclude
-                                         // us using non-primitive elements
-                                         // in components that we aren't
-                                         // interested in, however. make
-                                         // sure that all shape functions
-                                         // that are non-zero for the
-                                         // components we are interested in,
-                                         // are in fact primitive
-        for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
-          {
-            const std::vector<bool> &nonzero_component_array
-              = cell->get_fe().get_nonzero_components (i);
-            for (unsigned int c=0; c<n_components; ++c)
-              if ((nonzero_component_array[c] == true)
-                  &&
-                  (component_mask[c] == true))
-                Assert (cell->get_fe().is_primitive (i),
-                        ExcMessage ("This function can only deal with requested boundary "
-                                    "values that correspond to primitive (scalar) base "
-                                    "elements"));
-          }
-
-        typename MGDoFHandler<dim,spacedim>::face_iterator face = cell->face(face_no);
-        const unsigned char boundary_component = face->boundary_indicator();
-        if (function_map.find(boundary_component) != function_map.end())
-                                           // face is of the right component
-          {
-                                             // get indices, physical location and
-                                             // boundary values of dofs on this
-                                             // face
-            face_dofs.resize (fe.dofs_per_face);
-            face->get_mg_dof_indices (level, face_dofs);
-            if (fe_is_system)
-              {
-                                                 // enter those dofs
-                                                 // into the list that
-                                                 // match the
-                                                 // component
-                                                 // signature. avoid
-                                                 // the usual
-                                                 // complication that
-                                                 // we can't just use
-                                                 // *_system_to_component_index
-                                                 // for non-primitive
-                                                 // FEs
-                for (unsigned int i=0; i<face_dofs.size(); ++i)
-                  {
-                    unsigned int component;
-                    if (fe.is_primitive())
-                      component = fe.face_system_to_component_index(i).first;
-                    else
-                      {
-                                                         // non-primitive
-                                                         // case. make
-                                                         // sure that
-                                                         // this
-                                                         // particular
-                                                         // shape
-                                                         // function
-                                                         // _is_
-                                                         // primitive,
-                                                         // and get at
-                                                         // it's
-                                                         // component. use
-                                                         // usual
-                                                         // trick to
-                                                         // transfer
-                                                         // face dof
-                                                         // index to
-                                                         // cell dof
-
-                                                         // index
-                        const unsigned int cell_i
-                          = (dim == 1 ?
-                             i
-                             :
-                             (dim == 2 ?
-                              (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
-                              :
-                              (dim == 3 ?
-                               (i<4*fe.dofs_per_vertex ?
-                                i
-                                :
-                                (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
-                                 i+4*fe.dofs_per_vertex
-                                 :
-                                 i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
-                               :
-                               numbers::invalid_unsigned_int)));
-                        Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
-
-                                                         // make sure
-                                                         // that if
-                                                         // this is
-                                                         // not a
-                                                         // primitive
-                                                         // shape function,
-                                                         // then all
-                                                         // the
-                                                         // corresponding
-                                                         // components
-                                                         // in the
-                                                         // mask are
-                                                         // not set
+      
+                                      // field to store the indices
+      std::vector<unsigned int> local_dofs;
+      local_dofs.reserve (DoFTools::max_dofs_per_face(dof));
+      std::fill (local_dofs.begin (), local_dofs.end (),
+                DoFHandler<dim,spacedim>::invalid_dof_index);
+      
+      typename MGDoFHandler<dim,spacedim>::cell_iterator
+       cell = dof.begin(),
+       endc = dof.end();
+      for (; cell!=endc; ++cell)
+       for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+            ++face_no)
+         {
+           const FiniteElement<dim> &fe = cell->get_fe();
+           const unsigned int level = cell->level();
+           
+                                            // we can presently deal only with
+                                            // primitive elements for boundary
+                                            // values. this does not preclude
+                                            // us using non-primitive elements
+                                            // in components that we aren't
+                                            // interested in, however. make
+                                            // sure that all shape functions
+                                            // that are non-zero for the
+                                            // components we are interested in,
+                                            // are in fact primitive
+           for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+             {
+               const std::vector<bool> &nonzero_component_array
+                 = cell->get_fe().get_nonzero_components (i);
+               for (unsigned int c=0; c<n_components; ++c)
+                 if ((nonzero_component_array[c] == true)
+                     &&
+                     (component_mask[c] == true))
+                   Assert (cell->get_fe().is_primitive (i),
+                           ExcMessage ("This function can only deal with requested boundary "
+                                       "values that correspond to primitive (scalar) base "
+                                       "elements"));
+             }
+           
+           typename MGDoFHandler<dim,spacedim>::face_iterator face = cell->face(face_no);
+           const unsigned char boundary_component = face->boundary_indicator();
+           if (function_map.find(boundary_component) != function_map.end())
+                                              // face is of the right component
+             {
+                                                // get indices, physical location and
+                                                // boundary values of dofs on this
+                                                // face
+               local_dofs.resize (fe.dofs_per_face);
+               face->get_mg_dof_indices (level, local_dofs);
+               if (fe_is_system)
+                 {
+                                                    // enter those dofs
+                                                    // into the list that
+                                                    // match the
+                                                    // component
+                                                    // signature. avoid
+                                                    // the usual
+                                                    // complication that
+                                                    // we can't just use
+                                                    // *_system_to_component_index
+                                                    // for non-primitive
+                                                    // FEs
+                   for (unsigned int i=0; i<local_dofs.size(); ++i)
+                     {
+                       unsigned int component;
+                       if (fe.is_primitive())
+                         component = fe.face_system_to_component_index(i).first;
+                       else
+                         {
+                                                            // non-primitive
+                                                            // case. make
+                                                            // sure that
+                                                            // this
+                                                            // particular
+                                                            // shape
+                                                            // function
+                                                            // _is_
+                                                            // primitive,
+                                                            // and get at
+                                                            // it's
+                                                            // component. use
+                                                            // usual
+                                                            // trick to
+                                                            // transfer
+                                                            // face dof
+                                                            // index to
+                                                            // cell dof
+                           
+                                                            // index
+                           const unsigned int cell_i
+                             = (dim == 1 ?
+                                i
+                                :
+                                (dim == 2 ?
+                                 (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
+                                 :
+                                 (dim == 3 ?
+                                  (i<4*fe.dofs_per_vertex ?
+                                   i
+                                   :
+                                   (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
+                                    i+4*fe.dofs_per_vertex
+                                    :
+                                    i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
+                                  :
+                                  numbers::invalid_unsigned_int)));
+                           Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
+
+                                                            // make sure
+                                                            // that if
+                                                            // this is
+                                                            // not a
+                                                            // primitive
+                                                            // shape function,
+                                                            // then all
+                                                            // the
+                                                            // corresponding
+                                                            // components
+                                                            // in the
+                                                            // mask are
+                                                            // not set
 //                         if (!fe.is_primitive(cell_i))
 //                           for (unsigned int c=0; c<n_components; ++c)
 //                             if (fe.get_nonzero_components(cell_i)[c])
@@ -1279,22 +1318,23 @@ MGTools::make_boundary_list(
 // components. if shape function is non-primitive, then we will ignore
 // the result in the following anyway, otherwise there's only one
 // non-zero component which we will use
-                        component = (std::find (fe.get_nonzero_components(cell_i).begin(),
-                                                fe.get_nonzero_components(cell_i).end(),
-                                                true)
-                                     -
-                                     fe.get_nonzero_components(cell_i).begin());
-                      }
-
-                    if (component_mask[component] == true)
-                      boundary_indices[level].insert(face_dofs[i]);
-                  }
-              }
-            else
-             for (unsigned int i=0; i<face_dofs.size(); ++i)
-               boundary_indices[level].insert(face_dofs[i]);
-          }
-      }
+                           component = (std::find (fe.get_nonzero_components(cell_i).begin(),
+                                                   fe.get_nonzero_components(cell_i).end(),
+                                                   true)
+                                        -
+                                        fe.get_nonzero_components(cell_i).begin());
+                         }
+
+                       if (component_mask[component] == true)
+                         boundary_indices[level].insert(local_dofs[i]);
+                     }
+                 }
+               else
+                 for (unsigned int i=0; i<local_dofs.size(); ++i)
+                   boundary_indices[level].insert(local_dofs[i]);
+             }
+         }
+    }
 }
 
 #endif
index 558e32645e149b043d4e441ea7d499278563f704..c125ed4845440dabac506b4025554ac4392efd96 100644 (file)
@@ -3,7 +3,7 @@
 
 /*    $Id$       */
 /*                                                                */
-/*    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors */
+/*    Copyright (C) 2010 by the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */

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.