Fix a bug in function DoFCellAccessor::distribute_local_to_global()

doc/news/changes/minor/20240115tjhei

@@ -0,0 +1,7 @@
+Changed:
+parallel::distributed::Triangulation::prepare_coarsening_and_refinement()
+is now a collective call that exchanges coarsen/refinement flags on
+ghost cells. This makes adaptive refinement independent of the number
+of MPI ranks involved.
+<br>
+(Timo Heister, Quang Hoang, Vladimir Yushutin, 2024/01/15)

examples/step-25/step-25.cc

@@ -336,11 +336,12 @@ namespace Step25
     // First we assemble the Jacobian matrix $F'_h(U^{n,l})$, where $U^{n,l}$
     // is stored in the vector <code>solution</code> for convenience.
     system_matrix.copy_from(mass_matrix);
-    system_matrix.add(std::pow(time_step * theta, 2), laplace_matrix);
+    system_matrix.add(Utilities::fixed_power<2>(time_step * theta),
+                      laplace_matrix);
 
     SparseMatrix<double> tmp_matrix(sparsity_pattern);
     compute_nl_matrix(old_solution, solution, tmp_matrix);
-    system_matrix.add(std::pow(time_step * theta, 2), tmp_matrix);
+    system_matrix.add(Utilities::fixed_power<2>(time_step * theta), tmp_matrix);
 
     // Next we compute the right-hand side vector. This is just the
     // combination of matrix-vector products implied by the description of
@@ -351,17 +352,18 @@ namespace Step25
 
     mass_matrix.vmult(system_rhs, solution);
     laplace_matrix.vmult(tmp_vector, solution);
-    system_rhs.add(std::pow(time_step * theta, 2), tmp_vector);
+    system_rhs.add(Utilities::fixed_power<2>(time_step * theta), tmp_vector);
 
     mass_matrix.vmult(tmp_vector, old_solution);
     system_rhs.add(-1.0, tmp_vector);
     laplace_matrix.vmult(tmp_vector, old_solution);
-    system_rhs.add(std::pow(time_step, 2) * theta * (1 - theta), tmp_vector);
+    system_rhs.add(Utilities::fixed_power<2>(time_step) * theta * (1 - theta),
+                   tmp_vector);
 
     system_rhs.add(-time_step, M_x_velocity);
 
     compute_nl_term(old_solution, solution, tmp_vector);
-    system_rhs.add(std::pow(time_step, 2) * theta, tmp_vector);
+    system_rhs.add(Utilities::fixed_power<2>(time_step) * theta, tmp_vector);
 
     system_rhs *= -1.;
   }

examples/step-43/step-43.cc

@@ -300,7 +300,7 @@ namespace Step43
 
     const double numerator =
       2.0 * S * temp - S * S * (2.0 * S - 2.0 * viscosity * (1 - S));
-    const double denominator = std::pow(temp, 2.0);
+    const double denominator = Utilities::fixed_power<2>(temp);
 
     const double F_prime = numerator / denominator;
 

examples/step-44/step-44.cc

@@ -1862,7 +1862,7 @@ namespace Step44
             const double det_F_qp   = lqph[q_point]->get_det_F();
             const double J_tilde_qp = lqph[q_point]->get_J_tilde();
             const double the_error_qp_squared =
-              std::pow((det_F_qp - J_tilde_qp), 2);
+              Utilities::fixed_power<2>((det_F_qp - J_tilde_qp));
             const double JxW = fe_values.JxW(q_point);
 
             dil_L2_error += the_error_qp_squared * JxW;

examples/step-47/step-47.cc

@@ -130,7 +130,7 @@ namespace Step47
                            const unsigned int /*component*/ = 0) const override
 
       {
-        return 4 * std::pow(PI, 4.0) * std::sin(PI * p[0]) *
+        return 4 * Utilities::fixed_power<4>(PI) * std::sin(PI * p[0]) *
                std::sin(PI * p[1]);
       }
     };

examples/step-53/step-53.cc

@@ -263,9 +263,9 @@ namespace Step53
     const double th  = std::atan2(R * x(2), b * p);
     const double phi = std::atan2(x(1), x(0));
     const double theta =
-      std::atan2(x(2) + ep * ep * b * std::pow(std::sin(th), 3),
-                 (p -
-                  (ellipticity * ellipticity * R * std::pow(std::cos(th), 3))));
+      std::atan2(x(2) + ep * ep * b * Utilities::fixed_power<3>(std::sin(th)),
+                 (p - (ellipticity * ellipticity * R *
+                       Utilities::fixed_power<3>(std::cos(th)))));
     const double R_bar =
       R / (std::sqrt(1 - ellipticity * ellipticity * std::sin(theta) *
                            std::sin(theta)));

examples/step-55/step-55.cc

@@ -210,13 +210,15 @@ namespace Step55
                   std::exp(R_x * (-2 * std::sqrt(25.0 + 4 * pi2) + 10.0)) -
                 0.4 * pi2 * std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
                   std::cos(2 * R_y * pi) +
-                0.1 * std::pow(-std::sqrt(25.0 + 4 * pi2) + 5.0, 2) *
+                0.1 *
+                  Utilities::fixed_power<2>(-std::sqrt(25.0 + 4 * pi2) + 5.0) *
                   std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
                   std::cos(2 * R_y * pi);
     values[1] = 0.2 * pi * (-std::sqrt(25.0 + 4 * pi2) + 5.0) *
                   std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
                   std::sin(2 * R_y * pi) -
-                0.05 * std::pow(-std::sqrt(25.0 + 4 * pi2) + 5.0, 3) *
+                0.05 *
+                  Utilities::fixed_power<3>(-std::sqrt(25.0 + 4 * pi2) + 5.0) *
                   std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
                   std::sin(2 * R_y * pi) / pi;
 

examples/step-62/step-62.cc

@@ -402,10 +402,11 @@ namespace step62
             std::abs(p[1] - force_center[1]) < max_force_width_y / 2)
           {
             return max_force_amplitude *
-                   std::exp(-(std::pow(p[0] - force_center[0], 2) /
-                                (2 * std::pow(force_sigma_x, 2)) +
-                              std::pow(p[1] - force_center[1], 2) /
-                                (2 * std::pow(force_sigma_y, 2))));
+                   std::exp(
+                     -(Utilities::fixed_power<2>(p[0] - force_center[0]) /
+                         (2 * Utilities::fixed_power<2>(force_sigma_x)) +
+                       Utilities::fixed_power<2>(p[1] - force_center[1]) /
+                         (2 * Utilities::fixed_power<2>(force_sigma_y))));
           }
         else
           {
@@ -967,7 +968,7 @@ namespace step62
                   for (unsigned int j = 0; j < dofs_per_cell; ++j)
                     {
                       std::complex<double> matrix_sum = 0;
-                      matrix_sum += -std::pow(omega, 2) *
+                      matrix_sum += -Utilities::fixed_power<2>(omega) *
                                     quadrature_data.mass_coefficient[i][j];
                       matrix_sum += quadrature_data.stiffness_coefficient[i][j];
                       cell_matrix(i, j) += matrix_sum * quadrature_data.JxW;

examples/step-71/step-71.cc

@@ -2319,7 +2319,7 @@ namespace Step71
       // The first derivative of the saturation function, noting that
       // $\frac{d \tanh(x)}{dx} = \text{sech}^{2}(x)$.
       const double dtanh_two_h_dot_h_div_h_sat_squ =
-        std::pow(1.0 / std::cosh(two_h_dot_h_div_h_sat_squ), 2.0);
+        Utilities::fixed_power<2>(1.0 / std::cosh(two_h_dot_h_div_h_sat_squ));
       const Tensor<1, dim> dtwo_h_dot_h_div_h_sat_squ_dH =
         2.0 * 2.0 / (this->get_mu_e_h_sat() * this->get_mu_e_h_sat()) * H;
 

examples/step-85/step-85.cc

@@ -650,7 +650,7 @@ namespace Step85
                 const double      error_at_point =
                   solution_values.at(q) - analytical_solution.value(point);
                 error_L2_squared +=
-                  std::pow(error_at_point, 2) * fe_values->JxW(q);
+                  Utilities::fixed_power<2>(error_at_point) * fe_values->JxW(q);
               }
           }
       }

include/deal.II/base/tensor.h

@@ -29,6 +29,16 @@
 #  include <adolc/adouble.h> // Taped double
 #endif
 
+// boost::serialization::make_array used to be in array.hpp, but was
+// moved to a different file in BOOST 1.64
+#include <boost/version.hpp>
+#if BOOST_VERSION >= 106400
+#  include <boost/serialization/array_wrapper.hpp>
+#else
+#  include <boost/serialization/array.hpp>
+#endif
+
+
 #include <cmath>
 #include <ostream>
 
@@ -874,9 +884,7 @@ class Tensor
   /**
    * Array of tensors holding the subelements.
    */
-  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
-  // ... avoid a compiler warning in case of dim == 0 and ensure that the
-  // array always has positive size.
+  std::array<Tensor<rank_ - 1, dim, Number>, dim> values;
 
   /**
    * This constructor is for internal use. It provides a way
@@ -1289,7 +1297,7 @@ template <typename ArrayLike, std::size_t... indices>
 constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE
 Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
                                    std::index_sequence<indices...>)
-  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
+  : values{{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}}
 {
   static_assert(sizeof...(indices) == dim,
                 "dim should match the number of indices");
@@ -1356,7 +1364,9 @@ template <typename OtherNumber>
 constexpr DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number>::
 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
 {
-  return Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>(values);
+  Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> x;
+  std::copy(values.begin(), values.end(), x.values.begin());
+  return x;
 }
 
 
@@ -1813,7 +1823,10 @@ template <class Archive>
 inline void
 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
 {
-  ar &values;
+  if constexpr (rank_ > 1)
+    ar &values;
+  else
+    ar &boost::serialization::make_array(&values[0], dim);
 }
 
 

include/deal.II/distributed/tria.h

@@ -484,14 +484,17 @@ namespace parallel
        * Coarsen and refine the mesh according to refinement and coarsening
        * flags set.
        *
-       * Since the current processor only has control over those cells it owns
-       * (i.e. the ones for which <code>cell-@>subdomain_id() ==
-       * this-@>locally_owned_subdomain()</code>), refinement and coarsening
-       * flags are only respected for those locally owned cells. Flags may be
-       * set on other cells as well (and may often, in fact, if you call
-       * dealii::Triangulation::prepare_coarsening_and_refinement()) but will
-       * be largely ignored: the decision to refine the global mesh will only
-       * be affected by flags set on locally owned cells.
+       * Since the current processor only has control over those cells
+       * it owns (i.e. the ones for which <code>cell-@>subdomain_id()
+       * == this-@>locally_owned_subdomain()</code>), refinement and
+       * coarsening flags are only respected for those locally owned
+       * cells. Flags set on other cells will be ignored: the decision
+       * to refine the global mesh will only be affected by flags set
+       * on locally owned cells.
+       *
+       * This is a
+       * @ref GlossCollectiveOperation "collective operation"
+       * and needs to be called by all participating MPI ranks.
        *
        * @note This function by default partitions the mesh in such a way that
        * the number of cells on all processors is roughly equal. If you want
@@ -513,10 +516,16 @@ namespace parallel
       execute_coarsening_and_refinement() override;
 
       /**
-       * Override the implementation of prepare_coarsening_and_refinement from
-       * the base class. This is necessary if periodic boundaries are enabled
-       * and the level difference over vertices over the periodic boundary
-       * must not be more than 2:1.
+       * Prepare the triangulation for coarsening and refinement.
+       *
+       * This function performs necessary modifications of the
+       * coarsening and refinement flags to be consistent in parallel,
+       * to conform to smoothing flags set, and to conform to 2:1
+       * hanging node constraints.
+       *
+       * This is a
+       * @ref GlossCollectiveOperation "collective operation"
+       * and needs to be called by all participating MPI ranks.
        */
       virtual bool
       prepare_coarsening_and_refinement() override;
@@ -579,7 +588,9 @@ namespace parallel
       memory_consumption_p4est() const;
 
       /**
-       * A @ref GlossCollectiveOperation "collective operation" that produces a sequence of output files with
+       * A
+       * @ref GlossCollectiveOperation "collective operation"
+       * that produces a sequence of output files with
        * the given file base name that contain the mesh in VTK format.
        *
        * More than anything else, this function is useful for debugging the

include/deal.II/dofs/dof_accessor.templates.h

@@ -2957,7 +2957,7 @@ DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
     {
       global_matrix.add(dof_indices[i],
                         n_dofs,
-                        dof_indices.begin(),
+                        dof_indices.data(),
                         &local_matrix(i, 0));
       global_vector(dof_indices[i]) += local_vector(i);
     }

source/base/data_out_base.cc

@@ -4166,10 +4166,9 @@ namespace DataOutBase
                           h1(0) * h2(1) - h1(1) * h2(0);
 
                         // normalize Vector
-                        double norm =
-                          std::sqrt(std::pow(nrml[i * d1 + j * d2](0), 2.) +
-                                    std::pow(nrml[i * d1 + j * d2](1), 2.) +
-                                    std::pow(nrml[i * d1 + j * d2](2), 2.));
+                        double norm = std::hypot(nrml[i * d1 + j * d2](0),
+                                                 nrml[i * d1 + j * d2](1),
+                                                 nrml[i * d1 + j * d2](2));
 
                         if (nrml[i * d1 + j * d2](1) < 0)
                           norm *= -1.;