]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Skip hashes for comparing constraints in favor of a std::map with user-defined compar...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 17 Nov 2012 21:03:29 +0000 (21:03 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 17 Nov 2012 21:03:29 +0000 (21:03 +0000)
Fix error in dof_info.templates.h when hp elements are used but only one element is present.

Fix error in AVX detection.

git-svn-id: https://svn.dealii.org/trunk@27557 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/aclocal.m4
deal.II/configure
deal.II/include/deal.II/matrix_free/dof_info.templates.h
deal.II/include/deal.II/matrix_free/helper_functions.h
deal.II/include/deal.II/matrix_free/mapping_info.templates.h
deal.II/include/deal.II/matrix_free/matrix_free.h
deal.II/include/deal.II/matrix_free/matrix_free.templates.h

index 9fb7765dffad9fb53d4314e12a43986373f53d28..85da1c5535dfaab1924bbc9eea746ba63623f03c 100644 (file)
@@ -2210,7 +2210,7 @@ AC_DEFUN(DEAL_II_DETECT_VECTORIZATION_LEVEL, dnl
 [
   AC_LANG(C++)
   CXXFLAGS="$CXXFLAGSG"
-  dnl SSE2 check in debug mode
+  dnl SSE2 check
   AC_MSG_CHECKING(whether CPU supports SSE2)
   AC_TRY_RUN(
     [
@@ -2243,7 +2243,7 @@ AC_DEFUN(DEAL_II_DETECT_VECTORIZATION_LEVEL, dnl
     ],
     [
       AC_MSG_RESULT(yes)
-      dnl AVX check in debug mode
+      dnl AVX check
       AC_MSG_CHECKING(whether CPU supports AVX)
       AC_TRY_RUN(
       [
@@ -2271,6 +2271,7 @@ AC_DEFUN(DEAL_II_DETECT_VECTORIZATION_LEVEL, dnl
           if (ptr[i] != 5.0625)
             return_value = 1;
         _mm_free (data);
+        return return_value;
         }
       ],
       [
index 391072862f08563ecae302eeae4ce43990fefc54..5d843cb4e1414e40b63b2f2545d04fce97ac3c75 100755 (executable)
@@ -1,5 +1,5 @@
 #! /bin/sh
-# From configure.in Revision: 27122 .
+# From configure.in Revision: 27277 .
 # Guess values for system-dependent variables and create Makefiles.
 # Generated by GNU Autoconf 2.69 for deal.II 7.3.pre.
 #
@@ -814,7 +814,6 @@ enable_mpi
 enable_threads
 enable_shared
 enable_parser
-enable_mgcompatibility
 enable_compat_blocker
 with_boost
 with_petsc
@@ -1475,10 +1474,6 @@ Optional Features:
   --enable-parser         While switched on by default, this option allows to
                           switch off support for the function parser in the
                           contrib directory.
-  --enable-mgcompatibility
-                          Use preconditioner interface in MGSmootherRelaxation
-                          instead of the new interface using the function
-                          step. Defaults to disabled.
   --enable-compat-blocker=mapping
                           Block functions that implicitely assume a Q1 mapping
 
@@ -9859,17 +9854,6 @@ fi
 
 
 
-# Check whether --enable-mgcompatibility was given.
-if test "${enable_mgcompatibility+set}" = set; then :
-  enableval=$enable_mgcompatibility;  if test "x$enableval" = "xyes" ; then
-           { $as_echo "$as_me:${as_lineno-$LINENO}: result: enable multigrid compatibility mode" >&5
-$as_echo "enable multigrid compatibility mode" >&6; }
-
-$as_echo "#define DEAL_II_MULTIGRID_COMPATIBILITY 1" >>confdefs.h
-
-         fi
-fi
-
 
 
   # Check whether --enable-compat-blocker was given.
@@ -13583,6 +13567,7 @@ else
           if (ptr[i] != 5.0625)
             return_value = 1;
         _mm_free (data);
+        return return_value;
         }
 
 _ACEOF
index 8c680c96464649e855ad7cb311a990f57f745e97..6f2ee7ba7e729fbc444a672c6bdba1fd3ac8b2d8 100644 (file)
@@ -21,8 +21,17 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace internal
 {
-namespace MatrixFreeFunctions
-{
+  namespace MatrixFreeFunctions
+  {
+
+    struct ConstraintComparator
+    {
+      bool operator()(const std::pair<unsigned int,double> &p1,
+                      const std::pair<unsigned int,double> &p2) const
+      {
+        return p1.second < p2.second;
+      }
+    };
 
                                      /**
                                       * A struct that takes entries describing
@@ -52,135 +61,62 @@ namespace MatrixFreeFunctions
       unsigned short
       insert_entries (const std::vector<std::pair<unsigned int,double> > &entries);
 
-      std::vector<Number> constraint_pool_data;
-      std::vector<unsigned int> constraint_pool_row_index;
-      std::vector<std::pair<unsigned int, unsigned int> >  pool_locations;
-      std::vector<std::pair<Number,unsigned int> > constraint_entries;
+      std::vector<std::pair<unsigned int, double> > constraint_entries;
       std::vector<unsigned int> constraint_indices;
-      std::vector<double>      one_constraint;
-      HashValue                hashes;
+      std::pair<std::vector<Number>,unsigned int> next_constraint;
+      std::map<std::vector<Number>,unsigned int,FPArrayComparator<double> > constraints;
     };
 
 
     template <typename Number>
     ConstraintValues<Number>::ConstraintValues ()
       :
-      hashes (1.)
-    {
-      constraint_pool_row_index.push_back (0);
-    }
+      constraints(FPArrayComparator<Number>(1.))
+    {}
 
     template <typename Number>
     unsigned short
     ConstraintValues<Number>::
     insert_entries (const std::vector<std::pair<unsigned int,double> > &entries)
     {
-      unsigned int insert_position = deal_II_numbers::invalid_unsigned_int;
-
-      typedef std::vector<std::pair<unsigned int, unsigned int> >::iterator iter;
-      constraint_entries.resize(entries.size());
-      one_constraint.resize(entries.size());
-      constraint_indices.resize(entries.size());
-      for (unsigned int j=0;j<entries.size(); j++)
-        {
-          constraint_entries[j].first = entries[j].second;
-          constraint_entries[j].second = entries[j].first;
-        }
-
-      std::sort(constraint_entries.begin(),constraint_entries.end());
-      for (unsigned int j=0;j<entries.size(); j++)
+      next_constraint.first.resize(entries.size());
+      if (entries.size() > 0)
         {
+          constraint_indices.resize(entries.size());
+          constraint_entries = entries;
+          std::sort(constraint_entries.begin(), constraint_entries.end(),
+                    ConstraintComparator());
+          for (unsigned int j=0;j<constraint_entries.size(); j++)
+            {
                           // copy the indices of the constraint entries after
                           // sorting.
-          constraint_indices[j] = constraint_entries[j].second;
+              constraint_indices[j] = constraint_entries[j].first;
 
                           // one_constraint takes the weights of the
                           // constraint
-          one_constraint[j] = constraint_entries[j].first;
+              next_constraint.first[j] = constraint_entries[j].second;
+            }
         }
+      next_constraint.second = constraints.size();
 
                                 // check whether or not constraint is already
-                                // in pool
-      unsigned int hash_val = hashes(one_constraint);
-      std::pair<unsigned int,unsigned int> test (hash_val, 0);
-
-                                // Try to find a constraint in the pool with
-                                // the same hash value.
-      iter pos = std::lower_bound (pool_locations.begin(),
-                                   pool_locations.end(),
-                                   test);
-
-                                // If constraint has to be added, which will
-                                // be its no.
-      test.second = constraint_pool_row_index.size()-1;
-
-                                // Hash value larger than all the ones
-                                // before. We need to add it.
-      if (pos == pool_locations.end())
-        goto insert;
-
-                                // A constraint in the pool with the same hash
-                                // value identified.
-      else if (pos->first == test.first)
-        {
-          bool is_same = true;
-          while(is_same == true)
-            {
-              if(one_constraint.size()!=
-                 (constraint_pool_row_index[pos->second+1]-
-                  constraint_pool_row_index[pos->second]))
-                                // The constraints have different length, and
-                                // hence different.
-                is_same = false;
-              else
-                for (unsigned int q=0; q<one_constraint.size(); ++q)
-                                // check whether or not all weights are the
-                                // same.
-                  if (std::fabs(constraint_pool_data[constraint_pool_row_index
-                                                     [pos->second]+q]-
-                                one_constraint[q])>hashes.scaling)
-                    {
-                      is_same = false;
-                      break;
-                    }
-              if (is_same == false)
-                {
-                                // Try if there is another constraint with the
-                                // same hash value.
-                  ++pos;
-                  if (pos != pool_locations.end() && pos->first == test.first)
-                    is_same = true;
-                  else
-                    goto insert;
-                }
-              else
-                {
-                                // The constraint is the same as the
-                                // (pos->second)th in the pool. Add the
-                                // location of constraint in pool
-                  insert_position = pos->second;
-                  break;
-                }
-            }
-        }
+                                // in pool. the initial implementation
+                                // computed a hash value based on the truncated
+                                // array (to given accuracy around 1e-13) in
+                                // order to easily detect different arrays and
+                                // then made a fine-grained check when the
+                                // hash values were equal. this was quite
+                                // lenghty and now we use a std::map with a
+                                // user-defined comparator to compare floating
+                                // point arrays to a tolerance 1e-13.
+      std::pair<typename std::map<std::vector<Number>, unsigned int,
+                                  FPArrayComparator<Number> >::iterator,
+                bool> it = constraints.insert(next_constraint);
+      unsigned int insert_position = deal_II_numbers::invalid_unsigned_int;
+      if (it.second == false)
+        insert_position = it.first->second;
       else
-        {
-                                // A new constraint has been identified. It
-                                // needs to be added (at the position
-                                // test->second in pool_locations).
-        insert:
-          pool_locations.insert (pos, test);
-
-                                // Remember hash value and location of
-                                // constraint.
-          constraint_pool_data.insert (constraint_pool_data.end(),
-                                       one_constraint.begin(),
-                                       one_constraint.end());
-          constraint_pool_row_index.push_back (constraint_pool_data.size());
-
-                                // Add the location of constraint in pool.
-          insert_position = test.second;
-        }
+        insert_position = next_constraint.second;
 
                                 // we want to store the result as a short
                                 // variable, so we have to make sure that the
@@ -214,6 +150,9 @@ namespace MatrixFreeFunctions
     dofs_per_cell (dof_info_in.dofs_per_cell),
     dofs_per_face (dof_info_in.dofs_per_face),
     store_plain_indices (dof_info_in.store_plain_indices),
+    cell_active_fe_index (dof_info_in.cell_active_fe_index),
+    max_fe_index (dof_info_in.max_fe_index),
+    fe_index_conversion (dof_info_in.fe_index_conversion),
     ghost_dofs (dof_info_in.ghost_dofs)
   {}
 
@@ -229,9 +168,13 @@ namespace MatrixFreeFunctions
     ghost_dofs.clear();
     dofs_per_cell.clear();
     dofs_per_face.clear();
+    n_components = 0;
     row_starts_plain_indices.clear();
     plain_dof_indices.clear();
     store_plain_indices = false;
+    cell_active_fe_index.clear();
+    max_fe_index = 0;
+    fe_index_conversion.clear();
   }
 
 
@@ -251,7 +194,7 @@ namespace MatrixFreeFunctions
     const unsigned int n_owned     = last_owned - first_owned;
     std::pair<unsigned short,unsigned short> constraint_iterator (0,0);
 
-    unsigned int dofs_this_cell = (cell_active_fe_index.size() == 0) ?
+    unsigned int dofs_this_cell = (cell_active_fe_index.empty()) ?
       dofs_per_cell[0] : dofs_per_cell[cell_active_fe_index[cell_number]];
     for (unsigned int i=0; i<dofs_this_cell; i++)
       {
@@ -390,21 +333,6 @@ namespace MatrixFreeFunctions
               }
           }
       }
-
-#ifdef DEBUG
-    {
-      unsigned int n_dofs = 0;
-      const std::pair<unsigned short,unsigned short> * blb = begin_indicators(cell_number);
-      for (unsigned int j=0; j<row_length_indicators(cell_number); ++j)
-        {
-          n_dofs += blb[j].first;
-          n_dofs += constraint_values.constraint_pool_row_index[blb[j].second+1]
-            - constraint_values.constraint_pool_row_index[blb[j].second];
-        }
-      n_dofs += constraint_iterator.first;
-      AssertDimension(n_dofs, row_length_indices(cell_number));
-    }
-#endif
   }
 
 
@@ -535,7 +463,7 @@ namespace MatrixFreeFunctions
                                                  numbers::invalid_unsigned_int);
     const unsigned int n_boundary_cells = boundary_cells.size();
     for (unsigned int j=0; j<n_boundary_cells; ++j)
-      reverse_numbering[boundary_cells[j]] = 
+      reverse_numbering[boundary_cells[j]] =
         j + size_info.vectorization_length*size_info.boundary_cells_start;
     unsigned int counter = 0;
     unsigned int j = 0;
@@ -660,7 +588,7 @@ namespace MatrixFreeFunctions
     size_info.boundary_cells_end   = (size_info.boundary_cells_end -
                                       size_info.boundary_cells_start);
     size_info.boundary_cells_start = 0;
-    
+
     AssertDimension (counter, size_info.n_active_cells);
     renumbering = Utilities::invert_permutation (reverse_numbering);
   }
@@ -822,7 +750,7 @@ namespace MatrixFreeFunctions
     new_row_starts[size_info.n_macro_cells] =
       std_cxx1x::tuple<unsigned int,unsigned int,unsigned int>
       (new_dof_indices.size(), new_constraint_indicator.size(), 0);
-      
+
     AssertDimension(dof_indices.size(), new_dof_indices.size());
     AssertDimension(constraint_indicator.size(),
                     new_constraint_indicator.size());
@@ -984,7 +912,8 @@ namespace MatrixFreeFunctions
       for(counter=start_nonboundary*vectorization_length;
           counter<size_info.n_active_cells; counter++)
         {
-          renumbering_fe_index[cell_active_fe_index[renumbering[counter]]].
+          renumbering_fe_index[cell_active_fe_index.empty() ? 0 :
+                               cell_active_fe_index[renumbering[counter]]].
             push_back(renumbering[counter]);
         }
       counter = start_nonboundary * vectorization_length;
@@ -1010,14 +939,14 @@ namespace MatrixFreeFunctions
       if (task_info.block_size_last == 0)
         task_info.block_size_last = task_info.block_size;
     }
-    
+
                                     // assume that all FEs have the same
                                 // connectivity graph, so take the zeroth FE
     task_info.n_blocks = (size_info.n_macro_cells+task_info.block_size-1)/
       task_info.block_size;
     task_info.block_size_last = size_info.n_macro_cells-
       (task_info.block_size*(task_info.n_blocks-1));
-                                
+
                                     // create the connectivity graph with internal
                                 // blocking
     CompressedSimpleSparsityPattern connectivity;
@@ -1098,7 +1027,7 @@ namespace MatrixFreeFunctions
                         neighbor_neighbor_list.push_back(*neighbor);
                         partition_list[counter++] = *neighbor;
                       }
-                  }                  
+                  }
               }
             neighbor_list = neighbor_neighbor_list;
             neighbor_neighbor_list.resize(0);
@@ -1573,7 +1502,8 @@ namespace MatrixFreeFunctions
                         for(cell=counter-partition_counter; cell<counter; ++cell)
                           {
                             renumbering_fe_index
-                              [cell_active_fe_index[partition_partition_list
+                              [cell_active_fe_index.empty() ? 0 :
+                               cell_active_fe_index[partition_partition_list
                                                     [cell]]].
                               push_back(partition_partition_list[cell]);
                           }
@@ -1640,7 +1570,8 @@ namespace MatrixFreeFunctions
                               {
                                 unsigned int this_index = 0;
                                 if(hp_bool == true)
-                                  this_index = cell_active_fe_index[*neighbor];
+                                  this_index = cell_active_fe_index.empty() ? 0 :
+                                    cell_active_fe_index[*neighbor];
 
                                 // Only add this cell if we need more macro
                                 // cells in the current block or if there is a
@@ -2107,13 +2038,13 @@ namespace MatrixFreeFunctions
       (out, (row_starts.capacity()*sizeof(std_cxx1x::tuple<unsigned int,
                                           unsigned int, unsigned int>)));
     out << "       Memory dof indices:           ";
-    size_info.print_memory_statistics 
+    size_info.print_memory_statistics
       (out, MemoryConsumption::memory_consumption (dof_indices));
     out << "       Memory constraint indicators: ";
     size_info.print_memory_statistics
       (out, MemoryConsumption::memory_consumption (constraint_indicator));
     out << "       Memory plain indices:         ";
-    size_info.print_memory_statistics 
+    size_info.print_memory_statistics
       (out, MemoryConsumption::memory_consumption (row_starts_plain_indices)+
        MemoryConsumption::memory_consumption (plain_dof_indices));
     out << "       Memory vector partitioner:    ";
@@ -2169,7 +2100,7 @@ namespace MatrixFreeFunctions
   }
 
 
-} // end of namespace MatrixFreeFunctions
+  } // end of namespace MatrixFreeFunctions
 } // end of namespace internal
 
 DEAL_II_NAMESPACE_CLOSE
index 5fccb9d8b61f36a5875c6b32da30968042efa665..b9558fa3c0b038962f23c1db5a37353c2565f7d2 100644 (file)
@@ -22,8 +22,6 @@
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/vectorization.h>
 
-#include <boost/functional/hash.hpp>
-
 DEAL_II_NAMESPACE_OPEN
 
 
@@ -100,7 +98,7 @@ namespace MatrixFreeFunctions
                                  * to zero.
                                  */
     void clear();
-    
+
                                 /**
                                  * Prints minimum, average, and
                                  * maximal memory consumption over the
@@ -144,54 +142,38 @@ namespace MatrixFreeFunctions
 
                                 /**
                                  * Data type to identify cell type.
-                                 */ 
+                                 */
   enum CellType {cartesian=0, affine=1, general=2, undefined=3};
 
-  // ----------------- hash structure --------------------------------
 
-                                /**
-                                 * A class that is
-                                 * used to quickly find out whether two
-                                 * vectors of floating point numbers are the
-                                 * same without going through all the
-                                 * elements: store a hash value for each
-                                 * vector. Generate the
-                                 * hash value by a sum of all values
-                                 * multiplied by random numbers (cast to
-                                 * int). Of course, this is not a sure
-                                 * criterion and one must manually check for
-                                 * equality before this hash is telling
-                                 * something useful. However, inequalities are
-                                 * easily detected (unless roundoff spoils the
-                                 * hash function)
-                                 */
-  struct HashValue
+  /**
+   * A class that is used to compare floating point arrays (e.g. std::vectors,
+   * Tensor<1,dim>, etc.). The idea of this class is to consider two arrays as
+   * equal if they are the same within a given tolerance. We use this
+   * comparator class within an std::map<> of the given arrays. Note that this
+   * comparison operator does not satisfy all the mathematical properties one
+   * usually wants to have (consider e.g. the numbers a=0, b=0.1, c=0.2 with
+   * tolerance 0.15; the operator gives a<c, but neither of a<b? or b<c? is
+   * satisfied). This is not a problem in the use cases for this class, but be
+   * careful when using it in other contexts.
+   */
+  template<typename Number>
+  struct FPArrayComparator
   {
-                                // Constructor: sets the size of Number values
-                                // with the typical magnitude that is to be
-                                // expected.
-    HashValue (const double element_size = 1.);
-    
-                                // get hash value for a vector of floating
-                                // point numbers (which are assumed to be of
-                                // order of magnitude one). Do this by first
-                                // truncating everything that is smaller than
-                                // the scaling (in order to eliminate noise
-                                // from roundoff errors) and then calling the
-                                // boost hash function
-    unsigned int operator ()(const std::vector<double> &vec);
-
-                                // get hash value for a tensor of rank
-                                // two where the magnitude of the
-                                // entries is given by the parameter
-                                // weight
-    template <int dim, typename number>
-    unsigned int operator ()(const Tensor<2,dim,VectorizedArray<number> > 
-                             &input,
-                             const bool     is_diagonal);
-    
-    const double scaling;
+    FPArrayComparator (const Number scaling);
+
+    bool operator() (const std::vector<Number> &v1,
+                     const std::vector<Number> &v2) const;
+
+    template <int dim>
+    bool operator ()(const Tensor<1,dim,VectorizedArray<Number> > *t1,
+                     const Tensor<1,dim,VectorizedArray<Number> > *t2) const;
+
+    template <int dim>
+    bool operator ()(const Tensor<2,dim,VectorizedArray<Number> > *t1,
+                     const Tensor<2,dim,VectorizedArray<Number> > *t2) const;
+
+    Number tolerance;
   };
 
   // Note: Implementation in matrix_free.templates.h
index 9d1626f45a9e08af1301172e24433ece5a51d064..f36c058080eb54eddc6c6849e035a7975104ce18 100644 (file)
@@ -202,7 +202,7 @@ namespace MatrixFreeFunctions
               Utilities::fixed_power<dim>(n_q_points_1d[q]);
             current_data.n_q_points.push_back (n_q_points);
 
-            current_data.n_q_points_face.push_back 
+            current_data.n_q_points_face.push_back
               (Utilities::fixed_power<dim-1>(n_q_points_1d[q]));
             current_data.quadrature.push_back
               (Quadrature<dim>(quad[my_q][q]));
@@ -245,7 +245,7 @@ namespace MatrixFreeFunctions
                                 // considered together, this variable holds
                                 // the individual info of the last chunk of
                                 // cells
-        CellType cell_t [vectorization_length], 
+        CellType cell_t [vectorization_length],
           cell_t_prev [vectorization_length];
         for (unsigned int j=0; j<vectorization_length; ++j)
           cell_t_prev[j] = undefined;
@@ -287,11 +287,11 @@ namespace MatrixFreeFunctions
               current_data.jacobians_grad_upper.reserve (reserve_size);
           }
 
-                                // a hash structure that is used to detect
-                                // similarities between mapping data from one
-                                // cell to the next.
-        std::vector<std::pair<unsigned int, int> > hash_collection;
-        HashValue hash_value (jacobian_size);
+        FPArrayComparator<Number> comparator(jacobian_size);
+        std::map<const Tensor<1,dim,VectorizedArray<Number> >*, unsigned int,
+                 FPArrayComparator<Number> > cartesians(comparator);
+        std::map<const Tensor<2,dim,VectorizedArray<Number> >*, unsigned int,
+                 FPArrayComparator<Number> > affines(comparator);
 
                                 // loop over all cells
         for (unsigned int cell=0; cell<n_macro_cells; ++cell)
@@ -305,7 +305,7 @@ namespace MatrixFreeFunctions
               active_fe_index[cell] : 0;
             const unsigned int n_q_points = current_data.n_q_points[fe_index];
             if (fe_values[fe_index].get() == 0)
-              fe_values[fe_index].reset 
+              fe_values[fe_index].reset
                 (new FEValues<dim> (mapping, dummy_fe,
                                     current_data.quadrature[fe_index],
                                     update_flags_feval));
@@ -341,105 +341,73 @@ namespace MatrixFreeFunctions
                   if (cell_t[j] > most_general_type)
                     most_general_type = cell_t[j];
                 AssertIndexRange (most_general_type, 3);
+                unsigned int insert_position = numbers::invalid_unsigned_int;
 
                                 // Cartesian cell with diagonal Jacobian: only
                                 // insert the diagonal of the inverse and the
-                                // Jacobian determinant
-                unsigned int insert_position = numbers::invalid_unsigned_int;
-                typedef std::vector<std::pair<unsigned int,int> >::iterator iter;
+                                // Jacobian determinant. We do this by using
+                                // an std::map that collects pointers to all
+                                // Cartesian Jacobians. We need a pointer in
+                                // the std::map because it cannot store data
+                                // based on VectorizedArray (alignment
+                                // issue). We circumvent the problem by
+                                // temporarily filling the next value into the
+                                // cartesian_data field and, in case we did an
+                                // insertion, the data is already in the
+                                // correct place.
                 if (most_general_type == cartesian)
                   {
-                    std::pair<Tensor<1,dim,VectorizedArray<Number> >,
-                              VectorizedArray<Number> > new_entry;
+                    std::pair<const Tensor<1,dim,VectorizedArray<Number> >*,
+                              unsigned int> new_entry;
+                    new_entry.second = cartesian_data.size();
+                    Tensor<1,dim,VectorizedArray<Number> > cart;
                     for (unsigned int d=0; d<dim; ++d)
-                      new_entry.first[d] = data.const_jac[d][d];
-                    insert_position = cartesian_data.size();
-
-                                // check whether everything is the same as on
-                                // another cell before. find an insertion point
-                                // in the list of keys that we have
-                                // collected. put negative position so that
-                                // the insertion position from lower_bound is
-                                // that with the same key
-                    const unsigned int hash = hash_value.template operator()<dim,Number> (data.const_jac, true);
-                    std::pair<unsigned int,int> insertion (hash, -insert_position);
-                    iter pos = std::lower_bound (hash_collection.begin(),
-                                                 hash_collection.end(),
-                                                 insertion);
-
-                                // ok, found a data field with the same
-                                // key. check whether we really hit a
-                                // duplicate, i.e., whether the hash really
-                                // was effective
-                    bool duplicate = true;
-                    if (pos != hash_collection.end() &&
-                        pos->first == hash)
-                      {
-                        for (unsigned int d=0; d<dim; ++d)
-                          for (unsigned int j=0; j<vectorization_length; ++j)
-                            if (std::fabs(data.const_jac[d][d][j]-
-                                          cartesian_data[-pos->second].first[d][j])>
-                                hash_value.scaling)
-                              duplicate = false;
-                      }
-                    else
-                      duplicate = false;
-
-                                // if no duplicate, insert the data
-                    if (duplicate == false)
+                      cart[d] = data.const_jac[d][d];
+                    cartesian_data.push_back
+                      (std::pair<Tensor<1,dim,VectorizedArray<Number> >,
+                                 VectorizedArray<Number> >
+                       (cart, VectorizedArray<Number>()));
+                    new_entry.first = &cartesian_data[new_entry.second].first;
+
+                    std::pair<typename std::
+                              map<const Tensor<1,dim,VectorizedArray<Number> >*,
+                                  unsigned int,
+                                  FPArrayComparator<Number> >::iterator,
+                              bool> it = cartesians.insert(new_entry);
+                    if (it.second == false)
                       {
-                        hash_collection.insert (pos, insertion);
-                        cartesian_data.push_back (new_entry);
+                        insert_position = it.first->second;
+                        cartesian_data.resize(new_entry.second);
                       }
-                                // else, remember the position
                     else
-                      insert_position = -pos->second;
+                      insert_position = new_entry.second;
                   }
 
                                 // Constant Jacobian case. same strategy as
                                 // before, but with other data fields
                 else if (most_general_type == affine)
                   {
-                    insert_position = affine_data.size();
-
-                                // check whether everything is the same as on
-                                // the previous cell
-                    const unsigned int hash = hash_value.template operator()<dim,Number> (data.const_jac, false);
-                    std::pair<unsigned int,int> insertion (hash, -insert_position);
-                    iter pos = std::lower_bound (hash_collection.begin(),
-                                                 hash_collection.end(),
-                                                 insertion);
-
-                                // ok, found a data field with the same
-                                // key. check whether we really hit a
-                                // duplicate
-                    bool duplicate = true;
-                    if (pos != hash_collection.end() &&
-                        pos->first == hash)
+                    std::pair<const Tensor<2,dim,VectorizedArray<Number> >*,
+                              unsigned int> new_entry;
+                    new_entry.second = affine_data.size();
+                    affine_data.push_back
+                      (std::pair<Tensor<2,dim,VectorizedArray<Number> >,
+                                 VectorizedArray<Number> >
+                       (data.const_jac, VectorizedArray<Number>()));
+                    new_entry.first = &affine_data[new_entry.second].first;
+
+                    std::pair<typename std::
+                              map<const Tensor<2,dim,VectorizedArray<Number> >*,
+                                  unsigned int,
+                                  FPArrayComparator<Number> >::iterator,
+                              bool> it = affines.insert(new_entry);
+                    if (it.second == false)
                       {
-                        for (unsigned int d=0; d<dim; ++d)
-                          for (unsigned int e=0; e<dim; ++e)
-                            for (unsigned int j=0; j<vectorization_length; ++j)
-                              if (std::fabs(data.const_jac[d][e][j]-
-                                            affine_data[-pos->second].first[d]
-                                            [e][j])>
-                                  hash_value.scaling)
-                                duplicate = false;
+                        insert_position = it.first->second;
+                        affine_data.resize(new_entry.second);
                       }
                     else
-                      duplicate = false;
-
-                    if (duplicate == false)
-                      {
-                        hash_collection.insert (pos, insertion);
-                        affine_data.push_back 
-                          (std::pair<Tensor<2,dim,VectorizedArray<Number> >,
-                           VectorizedArray<Number> >(data.const_jac,
-                                                     make_vectorized_array 
-                                                     (Number(0.))));
-                      }
-                    else
-                      insert_position = -pos->second;
+                      insert_position = new_entry.second;
                   }
 
                                 // general cell case: first resize the data
@@ -598,7 +566,7 @@ namespace MatrixFreeFunctions
 
             if (update_flags & update_quadrature_points)
               {
-                                 // eventually we turn to the quadrature points
+                                // eventually we turn to the quadrature points
                                 // that we can compress in case we have
                                 // Cartesian cells. we also need to reorder
                                 // them into arrays of vectorized data types.
@@ -636,7 +604,7 @@ namespace MatrixFreeFunctions
                   }
               }
           } // end for ( cell < n_macro_cells )
-        current_data.rowstart_jacobians.push_back 
+        current_data.rowstart_jacobians.push_back
           (current_data.jacobians.size());
         current_data.rowstart_q_points[n_macro_cells] =
           current_data.quadrature_points.size();
@@ -711,7 +679,7 @@ namespace MatrixFreeFunctions
                                 // and we already have determined that this
                                 // cell is either Cartesian or with constant
                                 // Jacobian, we have nothing more to do.
-        if (my_q > 0 && (get_cell_type(cell) == cartesian 
+        if (my_q > 0 && (get_cell_type(cell) == cartesian
                          || get_cell_type(cell) == affine) )
           continue;
 
@@ -957,7 +925,7 @@ namespace MatrixFreeFunctions
     if (general_size_glob > 0)
       {
         out << "      Memory Jacobian data:          ";
-        size_info.print_memory_statistics 
+        size_info.print_memory_statistics
           (out, MemoryConsumption::memory_consumption (jacobians) +
            MemoryConsumption::memory_consumption (JxW_values));
         out << "      Memory second derivative data: ";
@@ -976,7 +944,7 @@ namespace MatrixFreeFunctions
     if (quad_size_glob > 0)
       {
         out << "      Memory quadrature points:      ";
-        size_info.print_memory_statistics 
+        size_info.print_memory_statistics
           (out, MemoryConsumption::memory_consumption (rowstart_q_points) +
            MemoryConsumption::memory_consumption (quadrature_points));
       }
@@ -990,10 +958,10 @@ namespace MatrixFreeFunctions
                                                          const SizeInfo &size_info) const
   {
     out << "    Cell types:                      ";
-    size_info.print_memory_statistics 
+    size_info.print_memory_statistics
       (out, MemoryConsumption::memory_consumption (cell_type));
     out << "    Memory transformations compr:    ";
-    size_info.print_memory_statistics 
+    size_info.print_memory_statistics
       (out, MemoryConsumption::memory_consumption (affine_data) +
        MemoryConsumption::memory_consumption (cartesian_data));
     for (unsigned int j=0; j<mapping_data_gen.size(); ++j)
index 266a77f58e4ce26518a1afc041849409e0f36883..bf75c4a7b8a6ed53160c08a671c5e69d79287737 100644 (file)
@@ -1011,6 +1011,12 @@ public:
   const internal::MatrixFreeFunctions::DoFInfo &
   get_dof_info (const unsigned int fe_component = 0) const;
 
+                                /**
+                                 * Returns the number of weights in the
+                                 * constraint pool.
+                                 */
+  unsigned int n_constraint_pool_entries() const;
+
                                 /**
                                  * Returns a pointer to the first
                                  * number in the constraint pool data
@@ -1330,6 +1336,16 @@ MatrixFree<dim,Number>::get_dof_info (unsigned int dof_index) const
 
 
 
+template <int dim, typename Number>
+inline
+unsigned int
+MatrixFree<dim,Number>::n_constraint_pool_entries() const
+{
+  return constraint_pool_row_index.size()-1;
+}
+
+
+
 template <int dim, typename Number>
 inline
 const Number*
index c26e375a722ebe01389a1895103669f8234a0c61..8faeee770d91f5311641eec99e91b67a312a03fc 100644 (file)
@@ -771,19 +771,38 @@ void MatrixFree<dim,Number>::initialize_indices
     AssertDimension (cell_level_index.size(),size_info.n_macro_cells*vectorization_length);
   }
 
-                                // set constraint pool and reorder the indices
-  constraint_pool_row_index =
-    constraint_values.constraint_pool_row_index;
-  constraint_pool_data.resize (constraint_values.constraint_pool_data.size());
-  std::copy (constraint_values.constraint_pool_data.begin(),
-             constraint_values.constraint_pool_data.end(),
-             constraint_pool_data.begin());
-  for (unsigned int no=0; no<n_fe; ++no)
+                                // set constraint pool from the std::map and
+                                // reorder the indices
+  typename std::map<std::vector<double>, unsigned int,
+           internal::MatrixFreeFunctions::FPArrayComparator<Number> >::iterator
+    it = constraint_values.constraints.begin(),
+    end = constraint_values.constraints.end();
+  std::vector<const std::vector<double>*>
+    constraints (constraint_values.constraints.size());
+  unsigned int length = 0;
+  for ( ; it != end; ++it)
+    {
+      AssertIndexRange(it->second, constraints.size());
+      constraints[it->second] = &it->first;
+      length += it->first.size();
+    }
+  constraint_pool_data.clear();
+  constraint_pool_data.reserve (length);
+  constraint_pool_row_index.reserve(constraint_values.constraints.size()+1);
+  constraint_pool_row_index.resize(1, 0);
+  for (unsigned int i=0; i<constraints.size(); ++i)
     {
-      dof_info[no].reorder_cells(size_info, renumbering,
-                                 constraint_pool_row_index,
-                                 irregular_cells, vectorization_length);
+      Assert(constraints[i] != 0, ExcInternalError());
+      constraint_pool_data.insert(constraint_pool_data.end(),
+                                  constraints[i]->begin(),
+                                  constraints[i]->end());
+      constraint_pool_row_index.push_back(constraint_pool_data.size());
     }
+  AssertDimension(constraint_pool_data.size(), length);
+  for (unsigned int no=0; no<n_fe; ++no)
+    dof_info[no].reorder_cells(size_info, renumbering,
+                               constraint_pool_row_index,
+                               irregular_cells, vectorization_length);
 
   indices_are_initialized = true;
 }
@@ -1070,56 +1089,69 @@ namespace internal
 
 
 
-  HashValue::HashValue (const double element_size)
-    :
-    scaling (element_size * std::numeric_limits<double>::epsilon() *
-             1024.)
-  {}
+    template <typename Number>
+    FPArrayComparator<Number>::FPArrayComparator (const Number scaling)
+      : 
+      tolerance (scaling * std::numeric_limits<double>::epsilon() * 1024.)
+    {}
 
 
 
-  unsigned int HashValue::operator ()(const std::vector<double> &vec)
-  {
-    std::vector<double> mod_vec(vec);
-    for (unsigned int i=0; i<mod_vec.size(); ++i)
-      mod_vec[i] -= fmod (mod_vec[i], scaling);
-    return static_cast<unsigned int>(boost::hash_range (mod_vec.begin(), mod_vec.end()));
-  }
+    template <typename Number>
+    bool
+    FPArrayComparator<Number>::operator() (const std::vector<Number> &v1,
+                                           const std::vector<Number> &v2) const
+    {
+      const unsigned int s1 = v1.size(), s2 = v2.size();
+      if (s1 < s2)
+        return true;
+      else if (s1 > s2)
+        return false;
+      else
+        for (unsigned int i=0; i<s1; ++i)
+          if (v1[i] < v2[i] - tolerance)
+            return true;
+          else if (v1[i] > v2[i] + tolerance)
+            return false;
+      return false;
+    }
 
 
-  template <int dim, typename number>
-  unsigned int HashValue::operator ()
-    (const Tensor<2,dim,VectorizedArray<number> > &input,
-     const bool     is_diagonal)
-  {
-    const unsigned int vectorization_length =
-      VectorizedArray<number>::n_array_elements;
 
-    if (is_diagonal)
-      {
-        number mod_tensor [dim][vectorization_length];
-        for (unsigned int i=0; i<dim; ++i)
-          for (unsigned int j=0; j<vectorization_length; ++j)
-            mod_tensor[i][j] = input[i][i][j] - fmod (input[i][i][j],
-                                                      number(scaling));
-        return static_cast<unsigned int>
-          (boost::hash_range(&mod_tensor[0][0],
-                             &mod_tensor[dim-1][vectorization_length-1]+1));
-      }
-    else
-      {
-        number mod_tensor [dim][dim][vectorization_length];
-        for (unsigned int i=0; i<dim; ++i)
-          for (unsigned int d=0; d<dim; ++d)
-            for (unsigned int j=0; j<vectorization_length; ++j)
-              mod_tensor[i][d][j] = input[i][d][j] - fmod (input[i][d][j],
-                                                           number(scaling));
-        return static_cast<unsigned int>(boost::hash_range
-                                         (&mod_tensor[0][0][0],
-                                          &mod_tensor[dim-1][dim-1][vectorization_length-1]+1));
-      }
-  }
+    template <typename Number>
+    template <int dim>
+    bool
+    FPArrayComparator<Number>::
+    operator ()(const Tensor<1,dim,VectorizedArray<Number> > *t1,
+                const Tensor<1,dim,VectorizedArray<Number> > *t2) const
+    {
+      for (unsigned int d=0; d<dim; ++d)
+        for (unsigned int k=0;k<VectorizedArray<Number>::n_array_elements;++k)
+          if ((*t1)[d][k] < (*t2)[d][k] - tolerance)
+            return true;
+          else if ((*t1)[d][k] > (*t2)[d][k] + tolerance)
+            return false;
+      return false;
+    }
+
+
 
+    template <typename Number>
+    template <int dim>
+    bool
+    FPArrayComparator<Number>::
+    operator ()(const Tensor<2,dim,VectorizedArray<Number> > *t1,
+                const Tensor<2,dim,VectorizedArray<Number> > *t2) const
+    {
+      for (unsigned int d=0; d<dim; ++d)
+        for (unsigned int e=0; e<dim; ++e)
+          for (unsigned int k=0;k<VectorizedArray<Number>::n_array_elements;++k)
+            if ((*t1)[d][e][k] < (*t2)[d][e][k] - tolerance)
+              return true;
+            else if ((*t1)[d][e][k] > (*t2)[d][e][k] + tolerance)
+              return false;
+      return false;
+    }
   }
 }
 

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.