[
AC_LANG(C++)
CXXFLAGS="$CXXFLAGSG"
- dnl SSE2 check in debug mode
+ dnl SSE2 check
AC_MSG_CHECKING(whether CPU supports SSE2)
AC_TRY_RUN(
[
],
[
AC_MSG_RESULT(yes)
- dnl AVX check in debug mode
+ dnl AVX check
AC_MSG_CHECKING(whether CPU supports AVX)
AC_TRY_RUN(
[
if (ptr[i] != 5.0625)
return_value = 1;
_mm_free (data);
+ return return_value;
}
],
[
#! /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.
#
enable_threads
enable_shared
enable_parser
-enable_mgcompatibility
enable_compat_blocker
with_boost
with_petsc
--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
-# 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.
if (ptr[i] != 5.0625)
return_value = 1;
_mm_free (data);
+ return return_value;
}
_ACEOF
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
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
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)
{}
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();
}
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++)
{
}
}
}
-
-#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
}
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;
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);
}
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());
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;
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;
neighbor_neighbor_list.push_back(*neighbor);
partition_list[counter++] = *neighbor;
}
- }
+ }
}
neighbor_list = neighbor_neighbor_list;
neighbor_neighbor_list.resize(0);
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]);
}
{
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
(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: ";
}
-} // end of namespace MatrixFreeFunctions
+ } // end of namespace MatrixFreeFunctions
} // end of namespace internal
DEAL_II_NAMESPACE_CLOSE
#include <deal.II/base/tensor.h>
#include <deal.II/base/vectorization.h>
-#include <boost/functional/hash.hpp>
-
DEAL_II_NAMESPACE_OPEN
* to zero.
*/
void clear();
-
+
/**
* Prints minimum, average, and
* maximal memory consumption over the
/**
* 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
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]));
// 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;
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)
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));
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
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.
}
}
} // 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();
// 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;
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: ";
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));
}
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)
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
+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*
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;
}
- 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;
+ }
}
}
projects / dealii.git / commitdiff