Merge branch 'dealii:master' into fixbug-DoFCellAccessor

cmake/setup_write_config.cmake

@@ -104,6 +104,8 @@ endif()
 if(NOT "${_instructions}" STREQUAL "")
   to_string(_string ${_instructions})
   _both(" (${_string})\n")
+else()
+  _both("\n")
 endif()
 
 if(CMAKE_CROSSCOMPILING)

doc/news/changes/minor/20230917SchreterMunch_2

@@ -1,4 +1,4 @@
-New: Add alternative interfaces to RemotePointEvaluation::evaluate_and_proecss and
+New: Add alternative interfaces to RemotePointEvaluation::evaluate_and_process and
 RemotePointEvaluation::process_and_evaluate.
 <br>
 (Magdalena Schreter, Peter Munch, 2023/09/17)

include/deal.II/base/tensor.h

@@ -1380,43 +1380,18 @@ Tensor<rank_, dim, Number>::Tensor(Tensor<rank_, dim, Number> &&other) noexcept
 }
 #  endif
 
-namespace internal
-{
-  namespace TensorSubscriptor
-  {
-    template <typename ArrayElementType, int dim>
-    constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE ArrayElementType &
-    subscript(ArrayElementType  *values,
-              const unsigned int i,
-              std::integral_constant<int, dim>)
-    {
-      AssertIndexRange(i, dim);
-      return values[i];
-    }
-
-    template <typename ArrayElementType>
-    constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE ArrayElementType &
-    subscript(ArrayElementType *dummy,
-              const unsigned int,
-              std::integral_constant<int, 0>)
-    {
-      Assert(
-        false,
-        ExcMessage(
-          "Cannot access elements of an object of type Tensor<rank,0,Number>."));
-      return *dummy;
-    }
-  } // namespace TensorSubscriptor
-} // namespace internal
 
 
 template <int rank_, int dim, typename Number>
 constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE
   typename Tensor<rank_, dim, Number>::value_type &
   Tensor<rank_, dim, Number>::operator[](const unsigned int i)
 {
-  return dealii::internal::TensorSubscriptor::subscript(
-    values, i, std::integral_constant<int, dim>());
+  Assert(dim != 0,
+         ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
+  AssertIndexRange(i, dim);
+
+  return values[i];
 }
 
 
@@ -1425,6 +1400,8 @@ constexpr DEAL_II_ALWAYS_INLINE
   DEAL_II_HOST_DEVICE const typename Tensor<rank_, dim, Number>::value_type &
   Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
 {
+  Assert(dim != 0,
+         ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
   AssertIndexRange(i, dim);
 
   return values[i];
@@ -1620,58 +1597,30 @@ constexpr inline DEAL_II_ALWAYS_INLINE
 }
 
 
-namespace internal
+
+template <int rank_, int dim, typename Number>
+template <typename OtherNumber>
+constexpr inline DEAL_II_ALWAYS_INLINE
+  DEAL_II_HOST_DEVICE Tensor<rank_, dim, Number> &
+  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
 {
-  namespace TensorImplementation
-  {
-    template <int rank,
-              int dim,
-              typename Number,
-              typename OtherNumber,
-              std::enable_if_t<
-                !std::is_integral<
-                  typename ProductType<Number, OtherNumber>::type>::value &&
-                  !std::is_same_v<Number, Differentiation::SD::Expression>,
-                int> = 0>
-    constexpr DEAL_II_HOST_DEVICE inline DEAL_II_ALWAYS_INLINE void
-    division_operator(Tensor<rank, dim, Number> (&t)[dim],
-                      const OtherNumber &factor)
+  if constexpr (std::is_integral<
+                  typename ProductType<Number, OtherNumber>::type>::value ||
+                std::is_same_v<Number, Differentiation::SD::Expression>)
     {
-      const Number inverse_factor = Number(1.) / factor;
       // recurse over the base objects
       for (unsigned int d = 0; d < dim; ++d)
-        t[d] *= inverse_factor;
+        values[d] /= s;
     }
-
-
-    template <int rank,
-              int dim,
-              typename Number,
-              typename OtherNumber,
-              std::enable_if_t<
-                std::is_integral<
-                  typename ProductType<Number, OtherNumber>::type>::value ||
-                  std::is_same_v<Number, Differentiation::SD::Expression>,
-                int> = 0>
-    constexpr DEAL_II_HOST_DEVICE inline DEAL_II_ALWAYS_INLINE void
-    division_operator(Tensor<rank, dim, Number> (&t)[dim],
-                      const OtherNumber &factor)
+  else
     {
-      // recurse over the base objects
+      // If we can, avoid division by multiplying by the inverse of the given
+      // factor:
+      const Number inverse_factor = Number(1.) / s;
       for (unsigned int d = 0; d < dim; ++d)
-        t[d] /= factor;
+        values[d] *= inverse_factor;
     }
-  } // namespace TensorImplementation
-} // namespace internal
 
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-constexpr inline DEAL_II_ALWAYS_INLINE
-  DEAL_II_HOST_DEVICE Tensor<rank_, dim, Number> &
-  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
-{
-  internal::TensorImplementation::division_operator(values, s);
   return *this;
 }
 
@@ -1723,12 +1672,31 @@ constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE
   typename numbers::NumberTraits<Number>::real_type
   Tensor<rank_, dim, Number>::norm_square() const
 {
-  typename numbers::NumberTraits<Number>::real_type s = internal::NumberType<
-    typename numbers::NumberTraits<Number>::real_type>::value(0.0);
-  for (unsigned int i = 0; i < dim; ++i)
-    s += values[i].norm_square();
-
-  return s;
+  if constexpr (dim == 0)
+    return internal::NumberType<
+      typename numbers::NumberTraits<Number>::real_type>::value(0.0);
+  else if constexpr (rank_ == 1)
+    {
+      // For rank-1 tensors, the square of the norm is simply the sum of
+      // squares of the elements:
+      typename numbers::NumberTraits<Number>::real_type s =
+        numbers::NumberTraits<Number>::abs_square(values[0]);
+      for (unsigned int i = 1; i < dim; ++i)
+        s += numbers::NumberTraits<Number>::abs_square(values[i]);
+
+      return s;
+    }
+  else
+    {
+      // For higher-rank tensors, the square of the norm is the sum
+      // of squares of sub-tensors
+      typename numbers::NumberTraits<Number>::real_type s =
+        values[0].norm_square();
+      for (unsigned int i = 1; i < dim; ++i)
+        s += values[i].norm_square();
+
+      return s;
+    }
 }
 
 
@@ -1789,45 +1757,6 @@ Tensor<rank_, dim, Number>::component_to_unrolled_index(
 
 
 
-namespace internal
-{
-  // unrolled_to_component_indices is instantiated from DataOut for dim==0
-  // and rank=2. Make sure we don't have compiler warnings.
-
-  template <int dim>
-  DEAL_II_HOST_DEVICE inline constexpr unsigned int
-  mod(const unsigned int x)
-  {
-    return x % dim;
-  }
-
-  template <>
-  DEAL_II_HOST_DEVICE inline unsigned int
-  mod<0>(const unsigned int x)
-  {
-    Assert(false, ExcInternalError());
-    return x;
-  }
-
-  template <int dim>
-  DEAL_II_HOST_DEVICE inline constexpr unsigned int
-  div(const unsigned int x)
-  {
-    return x / dim;
-  }
-
-  template <>
-  DEAL_II_HOST_DEVICE inline unsigned int
-  div<0>(const unsigned int x)
-  {
-    Assert(false, ExcInternalError());
-    return x;
-  }
-
-} // namespace internal
-
-
-
 template <int rank_, int dim, typename Number>
 constexpr inline TableIndices<rank_>
 Tensor<rank_, dim, Number>::unrolled_to_component_indices(const unsigned int i)
@@ -1837,17 +1766,28 @@ Tensor<rank_, dim, Number>::unrolled_to_component_indices(const unsigned int i)
   AssertIndexRange(i, dummy);
   (void)dummy;
 
-  TableIndices<rank_> indices;
-
-  unsigned int remainder = i;
-  for (int r = rank_ - 1; r >= 0; --r)
+  if constexpr (dim == 0)
     {
-      indices[r] = internal::mod<dim>(remainder);
-      remainder  = internal::div<dim>(remainder);
+      Assert(false,
+             ExcMessage(
+               "A tensor with dimension 0 does not store any elements. "
+               "There is no indexing that can address its elements."));
+      return {};
     }
-  Assert(remainder == 0, ExcInternalError());
+  else
+    {
+      TableIndices<rank_> indices;
 
-  return indices;
+      unsigned int remainder = i;
+      for (int r = rank_ - 1; r >= 0; --r)
+        {
+          indices[r] = remainder % dim;
+          remainder  = remainder / dim;
+        }
+      Assert(remainder == 0, ExcInternalError());
+
+      return indices;
+    }
 }
 
 

include/deal.II/matrix_free/fe_evaluation.h

@@ -8468,8 +8468,23 @@ FEFaceEvaluation<dim,
     {
       this->face_orientations[0] = 0;
       this->face_numbers[0]      = face_number;
-      for (unsigned int i = 0; i < n_lanes; ++i)
-        this->cell_ids[i] = cell_index * n_lanes + i;
+      if (this->matrix_free->n_active_entries_per_cell_batch(this->cell) ==
+          n_lanes)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (unsigned int i = 0; i < n_lanes; ++i)
+            this->cell_ids[i] = cell_index * n_lanes + i;
+        }
+      else
+        {
+          unsigned int i = 0;
+          for (; i <
+                 this->matrix_free->n_active_entries_per_cell_batch(this->cell);
+               ++i)
+            this->cell_ids[i] = cell_index * n_lanes + i;
+          for (; i < n_lanes; ++i)
+            this->cell_ids[i] = numbers::invalid_unsigned_int;
+        }
       for (unsigned int i = 0; i < n_lanes; ++i)
         this->face_ids[i] =
           this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,

include/deal.II/matrix_free/fe_point_evaluation.h

@@ -881,7 +881,8 @@ class FEPointEvaluation
    * integrated at the points.
    *
    * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
    *
    */
   template <std::size_t stride_view>
@@ -908,7 +909,8 @@ class FEPointEvaluation
    * integrated at the points.
    *
    * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
    *
    */
   void
@@ -938,7 +940,8 @@ class FEPointEvaluation
    * integrated at the points.
    *
    * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
    *
    */
   template <std::size_t stride_view>
@@ -970,7 +973,8 @@ class FEPointEvaluation
    * integrated at the points.
    *
    * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
    *
    */
   void

tests/examples/CMakeLists.txt.in

@@ -27,6 +27,12 @@ foreach(_diff_file ${_diff_files})
     set(_destination "${CMAKE_CURRENT_SOURCE_DIR}/${_destination}")
     file(CREATE_LINK ${_file} ${_destination} SYMBOLIC)
   endforeach()
+
+  # Create a symbolic link in the temporary source directory
+  # so that we can access files in the example folders:
+  file(CREATE_LINK "${DEAL_II_SOURCE_DIR}/examples/${_step}"
+    "${CMAKE_CURRENT_SOURCE_DIR}/${_step}" SYMBOLIC
+  )
 endforeach()
 
 #

tests/mpi/affine_constraints_get_view_01.cc

@@ -192,8 +192,15 @@ test()
     pressure_constraints.distribute(check_vector_2.block(1));
   }
 
-  // Now the assertion part: These vectors are the same:
-  Assert(check_vector_1 == check_vector_2, ExcInternalError());
+  // Finally check that these two vectors are the same. Note that
+  // when we compile with native optimizations we might have a slight
+  // difference in results:
+  {
+    TrilinosWrappers::MPI::BlockVector result = check_vector_1;
+    result -= check_vector_2;
+    Assert(result.l2_norm() / check_vector_1.l2_norm() < 1e-8,
+           ExcInternalError());
+  }
 
   deallog << "OK" << std::endl;
 }