From 886ad73b2b52bda871af61c5eaa8c7c59b2595b7 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Fri, 28 Dec 2012 08:16:20 +0000 Subject: [PATCH] The std::map that stored pointers to Jacobians could lead to a situation where the pointers became invalid. Store the actual data in the map now. git-svn-id: https://svn.dealii.org/trunk@27862 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/matrix_free/dof_info.templates.h | 1 + .../deal.II/matrix_free/helper_functions.h | 8 +- .../matrix_free/mapping_info.templates.h | 104 +++++++++--------- .../matrix_free/matrix_free.templates.h | 16 +-- 4 files changed, 67 insertions(+), 62 deletions(-) diff --git a/deal.II/include/deal.II/matrix_free/dof_info.templates.h b/deal.II/include/deal.II/matrix_free/dof_info.templates.h index 98c38e7b3f..a60093ecfa 100644 --- a/deal.II/include/deal.II/matrix_free/dof_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/dof_info.templates.h @@ -1002,6 +1002,7 @@ no_constraint: { // To start up, set the start_up cell to // partition and list all its neighbors. + AssertIndexRange(start_up, cell_partition.size()); cell_partition[start_up] = partition; neighbor_list.push_back(start_up); partition_list[counter++] = start_up; diff --git a/deal.II/include/deal.II/matrix_free/helper_functions.h b/deal.II/include/deal.II/matrix_free/helper_functions.h index 6011b6106b..630178b7b1 100644 --- a/deal.II/include/deal.II/matrix_free/helper_functions.h +++ b/deal.II/include/deal.II/matrix_free/helper_functions.h @@ -166,12 +166,12 @@ namespace internal const std::vector &v2) const; template - bool operator ()(const Tensor<1,dim,VectorizedArray > *t1, - const Tensor<1,dim,VectorizedArray > *t2) const; + bool operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, + const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const; template - bool operator ()(const Tensor<2,dim,VectorizedArray > *t1, - const Tensor<2,dim,VectorizedArray > *t2) const; + bool operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, + const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const; Number tolerance; }; diff --git a/deal.II/include/deal.II/matrix_free/mapping_info.templates.h b/deal.II/include/deal.II/matrix_free/mapping_info.templates.h index 2d9ab24941..be7a70d403 100644 --- a/deal.II/include/deal.II/matrix_free/mapping_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/mapping_info.templates.h @@ -287,11 +287,16 @@ end_set: current_data.jacobians_grad_upper.reserve (reserve_size); } + // we would like to put a Tensor<1,dim,VectorizedArray > as + // key into the std::map, but std::map allocation does not align the + // allocated memory correctly, so put it into a tensor of the + // correct length instead FPArrayComparator comparator(jacobian_size); - std::map > *, unsigned int, - FPArrayComparator > cartesians(comparator); - std::map > *, unsigned int, - FPArrayComparator > affines(comparator); + typedef Tensor<1,VectorizedArray::n_array_elements,Number> VEC_ARRAY; + std::map, unsigned int, + FPArrayComparator > cartesians(comparator); + std::map, unsigned int, + FPArrayComparator > affines(comparator); // loop over all cells for (unsigned int cell=0; cell > *, - unsigned int> new_entry; - new_entry.second = cartesian_data.size(); - Tensor<1,dim,VectorizedArray > cart; + std::pair,unsigned int> new_entry; + new_entry.second = cartesians.size(); for (unsigned int d=0; d >, - VectorizedArray > - (cart, VectorizedArray())); - new_entry.first = &cartesian_data[new_entry.second].first; - - std::pair > *, - unsigned int, - FPArrayComparator >::iterator, + for (unsigned int v=0; v::n_array_elements; ++v) + new_entry.first[d][v] = data.const_jac[d][d][v]; + + std::pair, + unsigned int, FPArrayComparator >::iterator, bool> it = cartesians.insert(new_entry); if (it.second == false) - { - insert_position = it.first->second; - cartesian_data.resize(new_entry.second); - } + insert_position = it.first->second; else insert_position = new_entry.second; } @@ -387,25 +381,18 @@ end_set: // before, but with other data fields else if (most_general_type == affine) { - std::pair > *, - unsigned int> new_entry; - new_entry.second = affine_data.size(); - affine_data.push_back - (std::pair >, - VectorizedArray > - (data.const_jac, VectorizedArray())); - new_entry.first = &affine_data[new_entry.second].first; - - std::pair > *, - unsigned int, - FPArrayComparator >::iterator, + std::pair,unsigned int> new_entry; + new_entry.second = affines.size(); + for (unsigned int d=0; d::n_array_elements; ++v) + new_entry.first[d][e][v] = data.const_jac[d][e][v]; + + std::pair, + unsigned int, FPArrayComparator >::iterator, bool> it = affines.insert(new_entry); if (it.second == false) - { - insert_position = it.first->second; - affine_data.resize(new_entry.second); - } + insert_position = it.first->second; else insert_position = new_entry.second; } @@ -615,20 +602,37 @@ end_set: // fields and compute the JxW value. if (my_q == 0) { - for (unsigned int i=0; i, + unsigned int>::iterator it = cartesians.begin(); + it != cartesians.end(); ++it) { - VectorizedArray det = cartesian_data[i].first[0]; - for (unsigned int d=1; d det = make_vectorized_array(1.); for (unsigned int d=0; d jac_d; + for (unsigned int v=0; + v::n_array_elements; ++v) + jac_d[v] = it->first[d][v]; + cartesian_data[it->second].first[d] = 1./jac_d; + det *= jac_d; + } + cartesian_data[it->second].second = std::abs(det); } - for (unsigned int i=0; i, + unsigned int>::iterator it = affines.begin(); + it != affines.end(); ++it) { - VectorizedArray det = determinant(affine_data[i].first); - affine_data[i].first = transpose(invert(affine_data[i].first)); - affine_data[i].second = std::abs(det); + Tensor<2,dim,VectorizedArray > jac; + for (unsigned int d=0; d::n_array_elements; ++v) + jac[d][e][v] = it->first[d][e][v]; + + affine_data[it->second].first = transpose(invert(jac)); + affine_data[it->second].second = std::abs(determinant(jac)); } } } diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h index be5e69d84f..c5e0d479e0 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h @@ -1121,14 +1121,14 @@ namespace internal template bool FPArrayComparator:: - operator ()(const Tensor<1,dim,VectorizedArray > *t1, - const Tensor<1,dim,VectorizedArray > *t2) const + operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, + const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const { for (unsigned int d=0; d::n_array_elements; ++k) - if ((*t1)[d][k] < (*t2)[d][k] - tolerance) + if ((t1)[d][k] < (t2)[d][k] - tolerance) return true; - else if ((*t1)[d][k] > (*t2)[d][k] + tolerance) + else if ((t1)[d][k] > (t2)[d][k] + tolerance) return false; return false; } @@ -1139,15 +1139,15 @@ namespace internal template bool FPArrayComparator:: - operator ()(const Tensor<2,dim,VectorizedArray > *t1, - const Tensor<2,dim,VectorizedArray > *t2) const + operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, + const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const { for (unsigned int d=0; d::n_array_elements; ++k) - if ((*t1)[d][e][k] < (*t2)[d][e][k] - tolerance) + if ((t1)[d][e][k] < (t2)[d][e][k] - tolerance) return true; - else if ((*t1)[d][e][k] > (*t2)[d][e][k] + tolerance) + else if ((t1)[d][e][k] > (t2)[d][e][k] + tolerance) return false; return false; } -- 2.39.5