]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
reverted some nonsensical changes on dof_tools.
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 15 May 2007 21:53:52 +0000 (21:53 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 15 May 2007 21:53:52 +0000 (21:53 +0000)
git-svn-id: https://svn.dealii.org/trunk@14673 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_tools.h
deal.II/deal.II/source/dofs/dof_tools.cc

index 672720d2daec6a63f6561d59c61bd7debf4aef96..c79a0f91e9228c6bb261aa346a97b2ee998feb27 100644 (file)
@@ -481,6 +481,46 @@ 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
index 7ebc61abb1d47858e593204d472cca111b6a4980..59c4b20e54e3b5a9e93605e20124f6220a07e5f9 100644 (file)
@@ -20,6 +20,7 @@
 #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>
@@ -168,6 +169,108 @@ 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>
@@ -964,7 +1067,9 @@ namespace internal
                ExcInternalError());
        Assert (fe2.dofs_per_line <= fe1.dofs_per_line,
                ExcInternalError());
-       Assert (fe2.dofs_per_quad <= fe1.dofs_per_quad,
+       Assert ((dim < 3)
+               ||
+               (fe2.dofs_per_quad <= fe1.dofs_per_quad),
                ExcInternalError());
 
                                         // the idea here is to designate as
@@ -1387,6 +1492,11 @@ namespace internal
                                         * 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
@@ -1426,10 +1536,36 @@ 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)
+                    if ((face_constraints(row,i) != 0)
+                       &&
+                       (std::fabs(face_constraints(row,i)) >= 1e-14*abs_sum))
                       {
 #ifdef WOLFGANG
                         std::cout << "   " << slave_dofs[row]
@@ -2489,17 +2625,58 @@ namespace internal
                                                       // (this is what
                                                       // happens in
                                                       // hp/hp_hanging_nodes_01
-                                                      // for example); check
-                                                      // the latter somewhat
-                                                      // crudely by comparing
-                                                      // fe names
-                     if (cell->get_fe().get_name() !=
-                         cell->neighbor(face)->get_fe().get_name())
+                                                      // 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))
                        {
-                         Assert (cell->get_fe().dofs_per_face == 0,
-                                 ExcNotImplemented());
-                         Assert (cell->neighbor(face)->get_fe().dofs_per_face == 0,
+                         Assert (cell->get_fe().dofs_per_face
+                                 ==
+                                 cell->neighbor(face)->get_fe().dofs_per_face,
                                  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;
@@ -3278,7 +3455,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
-  fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
+  std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
 
   const FiniteElement<dim> &fe   = dof_handler.get_fe();
 
@@ -4895,6 +5072,58 @@ 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>

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.