The std::map that stored pointers to Jacobians could lead to a situation where the...

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

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 98c38e7b3f7159fd796d9ac95f1f638aaa198b9d..a60093ecfa7e894b3e2b7457e2ea0df4fc057e89 100644 (file)
--- 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 6011b6106bab7e4dd357baaa592def110db61ebe..630178b7b1120618ba221abba1f51ff8957dd14f 100644 (file)
--- 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<Number> &v2) const;
 
       template <int dim>
-      bool operator ()(const Tensor<1,dim,VectorizedArray<Number> > *t1,
-                       const Tensor<1,dim,VectorizedArray<Number> > *t2) const;
+      bool operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+                       const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;
 
       template <int dim>
-      bool operator ()(const Tensor<2,dim,VectorizedArray<Number> > *t1,
-                       const Tensor<2,dim,VectorizedArray<Number> > *t2) const;
+      bool operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+                       const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::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 2d9ab249413bdb3210ad6a88cb2e559450e6662e..be7a70d40325e312a067edcf09c93432343cf49b 100644 (file)
--- 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<Number> > 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<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);
+          typedef Tensor<1,VectorizedArray<Number>::n_array_elements,Number> VEC_ARRAY;
+          std::map<Tensor<1,dim,VEC_ARRAY>, unsigned int,
+                   FPArrayComparator<Number> > cartesians(comparator);
+          std::map<Tensor<2,dim,VEC_ARRAY>, unsigned int,
+                   FPArrayComparator<Number> > affines(comparator);
 
           // loop over all cells
           for (unsigned int cell=0; cell<n_macro_cells; ++cell)
@@ -357,28 +362,17 @@ end_set:
                   // correct place.
                   if (most_general_type == cartesian)
                     {
-                      std::pair<const Tensor<1,dim,VectorizedArray<Number> > *,
-                          unsigned int> new_entry;
-                      new_entry.second = cartesian_data.size();
-                      Tensor<1,dim,VectorizedArray<Number> > cart;
+                      std::pair<Tensor<1,dim,VEC_ARRAY>,unsigned int> new_entry;
+                      new_entry.second = cartesians.size();
                       for (unsigned int d=0; d<dim; ++d)
-                        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,
+                        for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                          new_entry.first[d][v] = data.const_jac[d][d][v];
+
+                      std::pair<typename std::map<Tensor<1,dim,VEC_ARRAY>,
+                          unsigned int, FPArrayComparator<Number> >::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<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,
+                      std::pair<Tensor<2,dim,VEC_ARRAY>,unsigned int> new_entry;
+                      new_entry.second = affines.size();
+                      for (unsigned int d=0; d<dim; ++d)
+                        for (unsigned int e=0; e<dim; ++e)
+                          for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                            new_entry.first[d][e][v] = data.const_jac[d][e][v];
+
+                      std::pair<typename std::map<Tensor<2,dim,VEC_ARRAY>,
+                          unsigned int, FPArrayComparator<Number> >::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<cartesian_data.size(); ++i)
+              cartesian_data.resize(cartesians.size());
+              for (typename std::map<Tensor<1,dim,VEC_ARRAY>,
+                                     unsigned int>::iterator it = cartesians.begin();
+                   it != cartesians.end(); ++it)
                 {
-                  VectorizedArray<Number> det = cartesian_data[i].first[0];
-                  for (unsigned int d=1; d<dim; ++d)
-                    det *= cartesian_data[i].first[d];
+                  VectorizedArray<Number> det = make_vectorized_array<Number>(1.);
                   for (unsigned int d=0; d<dim; ++d)
-                    cartesian_data[i].first[d] = 1./cartesian_data[i].first[d];
-                  cartesian_data[i].second = std::abs(det);
+                    {
+                      VectorizedArray<Number> jac_d;
+                      for (unsigned int v=0;
+                           v<VectorizedArray<Number>::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<affine_data.size(); ++i)
+              affine_data.resize(affines.size());
+              for (typename std::map<Tensor<2,dim,VEC_ARRAY>,
+                                     unsigned int>::iterator it = affines.begin();
+                   it != affines.end(); ++it)
                 {
-                  VectorizedArray<Number> 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<Number> > jac;
+                  for (unsigned int d=0; d<dim; ++d)
+                    for (unsigned int e=0; e<dim; ++e)
+                      for (unsigned int v=0;
+                           v<VectorizedArray<Number>::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 be5e69d84fe74e1213202cc7ba70c7ce0b0ae206..c5e0d479e0a53b0d8973c5de940b63f1b553e2e8 100644 (file)
--- 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 <int dim>
     bool
     FPArrayComparator<Number>::
-    operator ()(const Tensor<1,dim,VectorizedArray<Number> > *t1,
-                const Tensor<1,dim,VectorizedArray<Number> > *t2) const
+    operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+                const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,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)
+          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 <int dim>
     bool
     FPArrayComparator<Number>::
-    operator ()(const Tensor<2,dim,VectorizedArray<Number> > *t1,
-                const Tensor<2,dim,VectorizedArray<Number> > *t2) const
+    operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+                const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,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)
+            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;
     }