]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fixed bug in find_cells_adjacent_to_vertex. Fixed behavior of count_dofs_per_block...
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 15 May 2007 16:05:01 +0000 (16:05 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 15 May 2007 16:05:01 +0000 (16:05 +0000)
git-svn-id: https://svn.dealii.org/trunk@14671 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_tools.h
deal.II/deal.II/source/dofs/dof_tools.cc
deal.II/deal.II/source/grid/grid_tools.cc
deal.II/doc/news/changes.html

index c1dfd56e6ade83834181fb5ac5edd6b996e0bfbe..672720d2daec6a63f6561d59c61bd7debf4aef96 100644 (file)
@@ -481,46 +481,6 @@ class DoFTools
     make_sparsity_pattern (const DH                              &dof,
                           const std::vector<std::vector<bool> > &mask,
                           SparsityPattern                       &sparsity_pattern);
-
-                                    /**
-                                     * Construct a sparsity pattern that
-                                     * allows coupling degrees of freedom on
-                                     * two different but related meshes.
-                                     *
-                                     * The idea is that if the two given
-                                     * DoFHandler objects correspond to two
-                                     * different meshes (and potentially to
-                                     * different finite elements used on
-                                     * these cells), but that if the two
-                                     * triangulations they are based on are
-                                     * derived from the same coarse mesh
-                                     * through hierarchical refinement, then
-                                     * one may set up a problem where one
-                                     * would like to test shape functions
-                                     * from one mesh against the shape
-                                     * functions from another mesh. In
-                                     * particular, this means that shape
-                                     * functions from a cell on the first
-                                     * mesh are tested against those on the
-                                     * second cell that are located on the
-                                     * corresponding cell; this
-                                     * correspondence is something that the
-                                     * IntergridMap class can determine.
-                                     *
-                                     * This function then constructs a
-                                     * sparsity pattern for which the degrees
-                                     * of freedom that represent the rows
-                                     * come from the first given DoFHandler,
-                                     * whereas the ones that correspond to
-                                     * columns come from the second
-                                     * DoFHandler.
-                                     */
-    template <class DH, class SparsityPattern>
-    static
-    void
-    make_sparsity_pattern (const DH        &dof_row,
-                          const DH        &dof_col,
-                          SparsityPattern &sparsity);
     
                                     /**
                                      * Create the sparsity pattern for
@@ -1103,8 +1063,17 @@ class DoFTools
                                      * velocities into one block, and
                                      * the pressure into another).
                                      *
-                                     * The result is returned in
-                                     * @p dofs_per_component.
+                                     * The result is returned in @p
+                                     * dofs_per_component. Note that
+                                     * the size of @p
+                                     * dofs_per_component needs to be
+                                     * enough to hold all the indices
+                                     * specified in @p
+                                     * target_component. If this is
+                                     * not the case, an assertion is
+                                     * thrown. The indices not
+                                     * targetted by target_components
+                                     * are left untouched.
                                      */
     template <int dim>
     static void
@@ -1123,7 +1092,11 @@ class DoFTools
                                      * counting is done by
                                      * blocks. See @ref GlossBlock
                                      * "blocks" in the glossary for
-                                     * details.
+                                     * details. Again the vectors are
+                                     * assumed to have the correct
+                                     * size before calling this
+                                     * function. If this is not the
+                                     * case, an assertion is thrown.
                                      */
     template <int dim>
     static void
@@ -1150,7 +1123,7 @@ class DoFTools
                                      * This function can be used when
                                      * different variables shall be
                                      * discretized on different
-                                     * grids, where one grid is
+                                     * grids, where one grid is
                                      * coarser than the other. This
                                      * idea might seem nonsensical at
                                      * first, but has reasonable
index a47c140e62706a711af58639ea6897c1bb64f748..f936d341246bac98d0827bedd7c11a3a5dd4946c 100644 (file)
@@ -20,7 +20,6 @@
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
 #include <grid/intergrid_map.h>
-#include <grid/grid_tools.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_constraints.h>
@@ -169,108 +168,6 @@ DoFTools::make_sparsity_pattern (
 
 
 
-template <class DH, class SparsityPattern>
-void
-DoFTools::make_sparsity_pattern (
-  const DH        &dof_row,
-  const DH        &dof_col,
-  SparsityPattern &sparsity)
-{
-  const unsigned int n_dofs_row = dof_row.n_dofs();
-  const unsigned int n_dofs_col = dof_col.n_dofs();
-
-  Assert (sparsity.n_rows() == n_dofs_row,
-          ExcDimensionMismatch (sparsity.n_rows(), n_dofs_row));
-  Assert (sparsity.n_cols() == n_dofs_col,
-          ExcDimensionMismatch (sparsity.n_cols(), n_dofs_col));
-
-
-  const std::list<std::pair<typename DH::cell_iterator,
-                            typename DH::cell_iterator> >
-    cell_list
-    = GridTools::get_finest_common_cells (dof_row, dof_col);
-
-
-  typename std::list<std::pair<typename DH::cell_iterator,
-                               typename DH::cell_iterator> >
-    ::const_iterator
-    cell_iter = cell_list.begin();
-
-  for (; cell_iter!=cell_list.end(); ++cell_iter)
-    {
-      const typename DH::cell_iterator cell_row = cell_iter->first;
-      const typename DH::cell_iterator cell_col = cell_iter->second;
-
-      if (!cell_row->has_children() && !cell_col->has_children())
-        {
-          const unsigned int dofs_per_cell_row =
-           cell_row->get_fe().dofs_per_cell;
-          const unsigned int dofs_per_cell_col =
-           cell_col->get_fe().dofs_per_cell;
-          std::vector<unsigned int>
-           local_dof_indices_row(dofs_per_cell_row);
-          std::vector<unsigned int>
-           local_dof_indices_col(dofs_per_cell_col);
-          cell_row->get_dof_indices (local_dof_indices_row);
-          cell_col->get_dof_indices (local_dof_indices_col);
-          for (unsigned int i=0; i<dofs_per_cell_row; ++i)
-            for (unsigned int j=0; j<dofs_per_cell_col; ++j)
-              sparsity.add (local_dof_indices_row[i],
-                            local_dof_indices_col[j]);
-        }
-      else if (cell_row->has_children())
-        {
-          const std::vector<typename DH::active_cell_iterator >
-           child_cells = GridTools::get_active_child_cells<DH> (cell_row);
-          for (unsigned int i=0; i<child_cells.size(); i++)
-            {
-             const typename DH::active_cell_iterator
-               cell_row_child = child_cells[i];
-              const unsigned int dofs_per_cell_row =
-               cell_row_child->get_fe().dofs_per_cell;
-              const unsigned int dofs_per_cell_col =
-               cell_col->get_fe().dofs_per_cell;
-              std::vector<unsigned int>
-               local_dof_indices_row(dofs_per_cell_row);
-              std::vector<unsigned int>
-               local_dof_indices_col(dofs_per_cell_col);
-              cell_row_child->get_dof_indices (local_dof_indices_row);
-              cell_col->get_dof_indices (local_dof_indices_col);
-              for (unsigned int i=0; i<dofs_per_cell_row; ++i)
-                for (unsigned int j=0; j<dofs_per_cell_col; ++j)
-                  sparsity.add (local_dof_indices_row[i],
-                                local_dof_indices_col[j]);
-            }
-        }
-      else
-        {
-          std::vector<typename DH::active_cell_iterator>
-           child_cells = GridTools::get_active_child_cells<DH> (cell_col);
-          for (unsigned int i=0; i<child_cells.size(); i++)
-            {
-             const typename DH::active_cell_iterator
-               cell_col_child = child_cells[i];
-              const unsigned int dofs_per_cell_row =
-               cell_row->get_fe().dofs_per_cell;
-              const unsigned int dofs_per_cell_col =
-               cell_col_child->get_fe().dofs_per_cell;
-              std::vector<unsigned int>
-               local_dof_indices_row(dofs_per_cell_row);
-              std::vector<unsigned int>
-               local_dof_indices_col(dofs_per_cell_col);
-              cell_row->get_dof_indices (local_dof_indices_row);
-              cell_col_child->get_dof_indices (local_dof_indices_col);
-              for (unsigned int i=0; i<dofs_per_cell_row; ++i)
-                for (unsigned int j=0; j<dofs_per_cell_col; ++j)
-                  sparsity.add (local_dof_indices_row[i],
-                                local_dof_indices_col[j]);
-            }
-        }
-    }
-}
-
-
-
 #if deal_II_dimension == 1
 
 template <class DH, class SparsityPattern>
@@ -1067,9 +964,7 @@ namespace internal
                ExcInternalError());
        Assert (fe2.dofs_per_line <= fe1.dofs_per_line,
                ExcInternalError());
-       Assert ((dim < 3)
-               ||
-               (fe2.dofs_per_quad <= fe1.dofs_per_quad),
+       Assert (fe2.dofs_per_quad <= fe1.dofs_per_quad,
                ExcInternalError());
 
                                         // the idea here is to designate as
@@ -1483,20 +1378,12 @@ namespace internal
 
       
                                        /**
-                                       * Copy constraints into a constraint
-                                       * matrix object.
-                                       *
-                                        * This function removes zero
+                                        * This method removes zero
                                         * constraints and those, which
                                         * constrain a DoF which was
                                         * already eliminated in one of
                                         * the previous steps of the hp
                                         * hanging node procedure.
-                                       *
-                                       * It also suppresses very small
-                                       * entries in the constraint matrix to
-                                       * avoid making the sparsity pattern
-                                       * fuller than necessary.
                                         */
 #ifdef DEAL_II_ANON_NAMESPACE_BUG
       static
@@ -1536,36 +1423,10 @@ namespace internal
                     }
 
               if (constraint_already_satisfied == false)
-                {
-                                                  // add up the absolute
-                                                  // values of all
-                                                  // constraints in this line
-                                                  // to get a measure of
-                                                  // their absolute size
-                 double abs_sum = 0;
-                 for (unsigned int i=0; i<n_master_dofs; ++i)
-                   abs_sum += std::abs (face_constraints(row,i));
-                 
-                                                  // then enter those
-                                                  // constraints that are
-                                                  // larger than
-                                                  // 1e-14*abs_sum. everything
-                                                  // else probably originated
-                                                  // from inexact inversion
-                                                  // of matrices and similar
-                                                  // effects. having those
-                                                  // constraints in here will
-                                                  // only lead to problems
-                                                  // because it makes
-                                                  // sparsity patterns fuller
-                                                  // than necessary without
-                                                  // producing any
-                                                  // significant effect
+                { 
                   constraints.add_line (slave_dofs[row]);
                   for (unsigned int i=0; i<n_master_dofs; ++i)
-                    if ((face_constraints(row,i) != 0)
-                       &&
-                       (std::fabs(face_constraints(row,i)) >= 1e-14*abs_sum))
+                    if (face_constraints(row,i) != 0)
                       {
 #ifdef WOLFGANG
                         std::cout << "   " << slave_dofs[row]
@@ -2625,58 +2486,17 @@ namespace internal
                                                       // (this is what
                                                       // happens in
                                                       // hp/hp_hanging_nodes_01
-                                                      // for example).
-                                                      //
-                                                      // another possibility
-                                                      // is what happens in
-                                                      // crash_13. there, we
-                                                      // have
-                                                      // FESystem(FE_Q(1),FE_DGQ(0))
-                                                      // vs. FESystem(FE_Q(1),FE_DGQ(1)).
-                                                      // neither of them
-                                                      // dominates the
-                                                      // other. the point is
-                                                      // that it doesn't
-                                                      // matter, since
-                                                      // hopefully for this
-                                                      // case, both sides'
-                                                      // dofs should have
-                                                      // been unified.
-                                                      //
-                                                      // make sure this is
-                                                      // actually true. this
-                                                      // actually only
-                                                      // matters, of course,
-                                                      // if either of the two
-                                                      // finite elements
-                                                      // actually do have
-                                                      // dofs on the face
-                     if ((cell->get_fe().dofs_per_face != 0)
-                         ||
-                         (cell->neighbor(face)->get_fe().dofs_per_face != 0))
+                                                      // for example); check
+                                                      // the latter somewhat
+                                                      // crudely by comparing
+                                                      // fe names
+                     if (cell->get_fe().get_name() !=
+                         cell->neighbor(face)->get_fe().get_name())
                        {
-                         Assert (cell->get_fe().dofs_per_face
-                                 ==
-                                 cell->neighbor(face)->get_fe().dofs_per_face,
+                         Assert (cell->get_fe().dofs_per_face == 0,
+                                 ExcNotImplemented());
+                         Assert (cell->neighbor(face)->get_fe().dofs_per_face == 0,
                                  ExcNotImplemented());
-
-                                                          // (ab)use the master
-                                                          // and slave dofs
-                                                          // arrays for a
-                                                          // moment here
-                         master_dofs.resize (cell->get_fe().dofs_per_face);
-                         cell->face(face)
-                           ->get_dof_indices (master_dofs,
-                                              cell->active_fe_index ());
-
-                         slave_dofs.resize (cell->neighbor(face)->get_fe().dofs_per_face);
-                         cell->face(face)
-                           ->get_dof_indices (slave_dofs,
-                                              cell->neighbor(face)->active_fe_index ());
-
-                         for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
-                           Assert (master_dofs[i] == slave_dofs[i],
-                                   ExcInternalError());
                        }
                      
                      break;
@@ -3455,7 +3275,7 @@ DoFTools::extract_hanging_node_dofs (const DoFHandler<2> &dof_handler,
   Assert(selected_dofs.size() == dof_handler.n_dofs(),
         ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
                                   // preset all values by false
-  std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
+  fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
 
   const FiniteElement<dim> &fe   = dof_handler.get_fe();
 
@@ -3676,9 +3496,6 @@ DoFTools::count_dofs_per_component (
 {
   const FiniteElement<dim>& fe = dof_handler.get_fe();
   const unsigned int n_components = fe.n_components();
-  dofs_per_component.resize (n_components);
-  std::fill (dofs_per_component.begin(), dofs_per_component.end(), 0U);
-  
                                   // If the empty vector was given as
                                   // default argument, set up this
                                   // vector as identity.
@@ -3691,6 +3508,21 @@ DoFTools::count_dofs_per_component (
   
   Assert(target_component.size()==n_components,
         ExcDimensionMismatch(target_component.size(),n_components));
+  
+
+                                  // Check that target components
+                                  // contains sensible information
+                                  // and reset the counters for
+                                  // dofs_per_component
+  unsigned int size_dst = dofs_per_component.size();
+  for(unsigned int i=0; i<n_components; ++i) {
+    Assert(target_component[i] < size_dst,
+          ExcIndexRange(target_component[i],
+                        0,
+                        size_dst));
+    dofs_per_component[target_component[i]] = 0;
+  }
+
 
                                   // special case for only one
                                   // component. treat this first
@@ -3698,11 +3530,11 @@ DoFTools::count_dofs_per_component (
                                   // computations
   if (n_components == 1)
     {
-      dofs_per_component[0] = dof_handler.n_dofs();
+      dofs_per_component[target_component[0]] = dof_handler.n_dofs();
       return;
     }
 
-      
+  
                                   // otherwise determine the number
                                   // of dofs in each component
                                   // separately. do so in parallel
@@ -3725,7 +3557,8 @@ DoFTools::count_dofs_per_component (
                                           dofs_in_component[i], false);
     };
   threads.join_all ();
-
+  
+  unsigned int total_count = 0;
                                   // next count what we got
   unsigned int component = 0;
   for (unsigned int b=0;b<fe.n_base_elements();++b)
@@ -3738,11 +3571,15 @@ DoFTools::count_dofs_per_component (
        {
          for (unsigned int dd=0;dd<d;++dd)
            {
-             if (base.is_primitive() || (!only_once || dd==0))
-             dofs_per_component[target_component[component]]
-               += std::count(dofs_in_component[component].begin(),
-                             dofs_in_component[component].end(),
-                             true);
+             if (base.is_primitive() || (!only_once || dd==0)) {
+               unsigned int this_count = 
+                 std::count(dofs_in_component[component].begin(),
+                            dofs_in_component[component].end(),
+                            true);
+               dofs_per_component[target_component[component]]
+                 += this_count;
+               total_count += this_count;
+             }
              ++component;
            }
        }
@@ -3752,12 +3589,11 @@ DoFTools::count_dofs_per_component (
                                   // only valid if the finite element
                                   // is actually primitive, so
                                   // exclude other elements from this
+
   Assert (!dof_handler.get_fe().is_primitive()
           ||
-          (std::accumulate (dofs_per_component.begin(),
-                            dofs_per_component.end(), 0U)
-           == dof_handler.n_dofs()),
-         ExcInternalError());
+          (total_count == dof_handler.n_dofs()),
+         ExcDimensionMismatch(total_count, dof_handler.n_dofs()));
 }
 
 
@@ -3770,8 +3606,7 @@ DoFTools::count_dofs_per_block (
 {
   const FiniteElement<dim>& fe = dof_handler.get_fe();
   const unsigned int n_blocks = fe.n_blocks();
-  dofs_per_block.resize (n_blocks);
-  std::fill (dofs_per_block.begin(), dofs_per_block.end(), 0U);
+  unsigned int size_dst = dofs_per_block.size();
   
                                   // If the empty vector was given as
                                   // default argument, set up this
@@ -3785,14 +3620,22 @@ DoFTools::count_dofs_per_block (
   
   Assert(target_block.size()==n_blocks,
         ExcDimensionMismatch(target_block.size(),n_blocks));
-
+  
+                                  // Check constistency and reset
+                                  // counters
+  for(unsigned int i=0; i<n_blocks; ++i) {
+    Assert(target_block[i] < size_dst,
+          ExcIndexRange(target_block[i], 0, size_dst));
+    dofs_per_block[target_block[i]] = 0;
+  }
+  
                                   // special case for only one
                                   // block. treat this first
                                   // since it does not require any
                                   // computations
   if (n_blocks == 1)
     {
-      dofs_per_block[0] = dof_handler.n_dofs();
+      dofs_per_block[target_block[0]] = dof_handler.n_dofs();
       return;
     } 
                                   // otherwise determine the number
@@ -5049,58 +4892,6 @@ DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
  const Table<2,Coupling>&,
  CompressedBlockSparsityPattern&);
 
-
-template void
-DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
-                               SparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof_row,
- const DoFHandler<deal_II_dimension> &dof_col,
- SparsityPattern    &sparsity);
-template void
-DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
-                               CompressedSparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof_row,
- const DoFHandler<deal_II_dimension> &dof_col,
- CompressedSparsityPattern    &sparsity);
-template void
-DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
-                               BlockSparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof_row,
- const DoFHandler<deal_II_dimension> &dof_col,
- BlockSparsityPattern                &sparsity);
-template void
-DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
-                               CompressedBlockSparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof_row,
- const DoFHandler<deal_II_dimension> &dof_col,
- CompressedBlockSparsityPattern      &sparsity);
-
-template void
-DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
-                               SparsityPattern>
-(const hp::DoFHandler<deal_II_dimension> &dof_row,
- const hp::DoFHandler<deal_II_dimension> &dof_col,
- SparsityPattern    &sparsity);
-template void
-DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
-                               CompressedSparsityPattern>
-(const hp::DoFHandler<deal_II_dimension> &dof_row,
- const hp::DoFHandler<deal_II_dimension> &dof_col,
- CompressedSparsityPattern    &sparsity);
-template void
-DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
-                               BlockSparsityPattern>
-(const hp::DoFHandler<deal_II_dimension> &dof_row,
- const hp::DoFHandler<deal_II_dimension> &dof_col,
- BlockSparsityPattern                &sparsity);
-template void
-DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
-                               CompressedBlockSparsityPattern>
-(const hp::DoFHandler<deal_II_dimension> &dof_row,
- const hp::DoFHandler<deal_II_dimension> &dof_col,
- CompressedBlockSparsityPattern      &sparsity);
-
-
 // #if deal_II_dimension > 1
 template void
 DoFTools::make_boundary_sparsity_pattern<DoFHandler<deal_II_dimension>,SparsityPattern>
index c30d15faa2182d0319bc23c730b964c4f28e7326..6c0bec81b51f26e36e900b966a2dce00419f7cd4 100644 (file)
@@ -625,7 +625,7 @@ GridTools::find_cells_adjacent_to_vertex(const Container<dim> &container,
                                    // all faces to which this vertex
                                    // belongs, check the level of
                                    // the neighbor, and if it is coarser,
-                                   // then check whether the vortex is
+                                   // then check whether the vertex is
                                    // part of that neighbor or not.
                for (unsigned vface = 0; vface < dim; vface++)
                   {
@@ -661,7 +661,7 @@ GridTools::find_cells_adjacent_to_vertex(const Container<dim> &container,
                             {
                               bool found = false;
                               for (unsigned v=0; v<GeometryInfo<dim>::vertices_per_cell; v++)
-                                 if (cell->vertex_index(v) == (int)vertex)
+                                 if (nb->vertex_index(v) == (int)vertex)
                                   {
                                     found = true;
                                     break;
index eee2283bb3fc3b056027028707a9b72cbfaf6cb9..d8cf807d61fcf37c54dc4558f0d140e1a1c0b8c1 100644 (file)
@@ -1033,6 +1033,23 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+   <li> <p>Fixed: the function
+       <code>GridTools::find_cells_adjacent_to_vertex</code> 
+       was not detecting properly the coarse cells adjacent to
+       refined cells.
+       <br>
+       (Luca Heltai 2007/05/15)
+       </p>
+
+   <li> <p>Fixed: the two tools
+       <code>DoFTools::count_dofs_per_component</code> and
+       <code>DoFTools::count_dofs_per_block</code> where changing the
+       size of the destination vector. Consistently with (most of) the 
+       rest of the library, now the vectors are expected to be the 
+       right size before calling these functions. 
+       <br>
+       (Luca Heltai 2007/05/15)
+       </p>
 
   <li> <p>New: The classes <code
   class="class">MGTransferBlockSelect</code> and <code

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.