Merge pull request adamantine-sim#479 from Rombur/continuous_five_axis

Add support for curved scan path when using five-axis

Author: Rombur

source/ElectronBeamHeatSource.cc

@@ -28,7 +28,7 @@ void ElectronBeamHeatSource<dim>::update_time(double time)
       this->_scan_path.get_power_modifier(time);
   this->_source_on = (segment_power_modifier > 0.0);
   dealii::Point<3> const &path = this->_scan_path.value(time);
-  _quaternion = this->_scan_path.get_current_quaternion();
+  _quaternion = this->_scan_path.get_quaternion(time);
   // Copy the scalar value of path to the vectorized data type
   for (unsigned int d = 0; d < 3; ++d)
   {

source/GoldakHeatSource.cc

@@ -28,7 +28,7 @@ void GoldakHeatSource<dim>::update_time(double time)
       this->_scan_path.get_power_modifier(time);
   this->_source_on = (segment_power_modifier > 0.0);
   dealii::Point<3> const &path = this->_scan_path.value(time);
-  _quaternion = this->_scan_path.get_current_quaternion();
+  _quaternion = this->_scan_path.get_quaternion(time);
   // Copy the scalar value of path to the vectorized data type
   for (unsigned int d = 0; d < 3; ++d)
   {

source/Quaternion.cc

@@ -8,6 +8,7 @@
 #include <deal.II/base/signaling_nan.h>
 
 #include <cmath>
+#include <numeric>
 
 namespace adamantine
 {
@@ -81,6 +82,18 @@ Quaternion &Quaternion::operator*=(Quaternion const &other)
   return *this;
 }
 
+Quaternion &Quaternion::operator*=(double const &scalar)
+{
+  for (int i = 0; i < 4; ++i)
+  {
+    _quaternion[i] *= scalar;
+  }
+
+  build_rotation_matrices();
+
+  return *this;
+}
+
 Quaternion &Quaternion::operator/=(Quaternion const &other)
 {
   ASSERT(_is_valid, "");
@@ -162,4 +175,10 @@ void Quaternion::build_rotation_matrices()
     }
   }
 }
+
+double dot_product(Quaternion const &p, Quaternion const &q)
+{
+  return std::inner_product(p._quaternion.begin(), p._quaternion.end(),
+                            q._quaternion.begin(), 0.);
+}
 } // namespace adamantine

source/Quaternion.hh

@@ -1,4 +1,4 @@
-/* SPDX-FileCopyrightText: Copyright (c) 2025, the adamantine authors.
+/* SPDX-FileCopyrightText: Copyright (c) 2025 - 2026, the adamantine authors.
  * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
  */
 
@@ -77,6 +77,11 @@ public:
    */
   Quaternion &operator*=(Quaternion const &other);
 
+  /**
+   * Multiply the quaternion with a scalar.
+   */
+  Quaternion &operator*=(double const &scalar);
+
   /**
    * Multiply the quaternions using the inverse of @p other.
    */
@@ -88,6 +93,11 @@ public:
    */
   void pow(double const exp);
 
+  /**
+   * Compute the dot product of two quaternions.
+   */
+  friend double dot_product(Quaternion const &p, Quaternion const &q);
+
 private:
   /**
    * Build the different rotation matrices.

source/ScanPath.cc

@@ -196,24 +196,27 @@ void ScanPath::load_event_series_scan_path()
 
 void ScanPath::update_current_segment_info(
     double time, bool save_segment, unsigned int &current_segment,
-    dealii::Point<3> &segment_start_point, double &segment_start_time) const
+    dealii::Point<3> &segment_start_point, double &segment_start_time,
+    Quaternion &segment_start_rotation) const
 {
   // Get to the correct segment
   current_segment = _old_segment;
   while (time > _segment_list[current_segment].end_time)
   {
     ++current_segment;
   }
-  // Update the start position and time for the current segment
+  // Update the start position, time, and rotation for the current segment
   if (current_segment > 0)
   {
     segment_start_time = _segment_list[current_segment - 1].end_time;
     segment_start_point = _segment_list[current_segment - 1].end_point;
+    segment_start_rotation = _segment_list[current_segment - 1].end_rotation;
   }
   else
   {
     segment_start_time = 0.0;
     segment_start_point = _segment_list[current_segment].end_point;
+    segment_start_rotation = _segment_list[current_segment].end_rotation;
   }
 
   if (save_segment)
@@ -237,9 +240,11 @@ dealii::Point<3> ScanPath::value(double const time, bool save_segment) const
   // Get to the correct segment
   dealii::Point<3> segment_start_point;
   double segment_start_time = 0.0;
+  Quaternion segment_start_rotation;
   unsigned int current_segment = 0;
   update_current_segment_info(time, save_segment, current_segment,
-                              segment_start_point, segment_start_time);
+                              segment_start_point, segment_start_time,
+                              segment_start_rotation);
 
   // Calculate the position in the direction given by "component"
   dealii::Point<3> position =
@@ -261,9 +266,10 @@ double ScanPath::get_power_modifier(double const time) const
   // Get to the correct segment
   dealii::Point<3> segment_start_point;
   double segment_start_time = 0.0;
+  Quaternion segment_start_rotation;
   unsigned int current_segment = 0;
   update_current_segment_info(time, true, current_segment, segment_start_point,
-                              segment_start_time);
+                              segment_start_time, segment_start_rotation);
 
   return _segment_list[current_segment].power_modifier;
 }
@@ -278,19 +284,61 @@ dealii::Point<3> ScanPath::rotate(double const time,
     return point;
   }
 
+  return get_quaternion(time).rotate(point);
+}
+
+Quaternion ScanPath::get_quaternion(double const time) const
+{
+  if (!_five_axis)
+  {
+    return Quaternion(1.0, 0.0, 0.0, 0.0);
+  }
+
+  // If the current time is after the scan path data is over, return the
+  // end rotation.
+  if (time > _segment_list.back().end_time)
+  {
+    return _segment_list.back().end_rotation;
+  }
+
   // Get to the correct segment
   dealii::Point<3> segment_start_point;
   double segment_start_time = 0.0;
+  Quaternion segment_start_rotation;
   unsigned int current_segment = 0;
   update_current_segment_info(time, false, current_segment, segment_start_point,
-                              segment_start_time);
+                              segment_start_time, segment_start_rotation);
 
-  return _segment_list[current_segment].end_rotation.rotate(point);
-}
+  // If the start and end rotation are the same, return the start rotation
+  if (segment_start_rotation == _segment_list[current_segment].end_rotation)
+  {
+    return segment_start_rotation;
+  }
 
-Quaternion ScanPath::get_current_quaternion() const
-{
-  return _segment_list[_old_segment].end_rotation;
+  // Compute the interpolated rotation
+  double const duration =
+      _segment_list[current_segment].end_time - segment_start_time;
+  if (duration < 1e-12)
+  {
+    return _segment_list[current_segment].end_rotation;
+  }
+
+  Quaternion interpolated_rotation =
+      _segment_list[current_segment].end_rotation;
+
+  // If the dot product of the initial quaternion and the final rotation is
+  // negative, slerp will take the long way instead of the short path. In that
+  // case, we need to negate one of the quaternion.
+  if (dot_product(segment_start_rotation, interpolated_rotation) < 0.)
+  {
+    interpolated_rotation *= -1.;
+  }
+
+  interpolated_rotation /= segment_start_rotation;
+  interpolated_rotation.pow((time - segment_start_time) / duration);
+  interpolated_rotation *= segment_start_rotation;
+
+  return interpolated_rotation;
 }
 
 std::vector<ScanPathSegment> ScanPath::get_segment_list() const

source/ScanPath.hh

@@ -97,9 +97,10 @@ public:
                           dealii::Point<3> const &point) const;
 
   /**
-   * Return the quaternion associated with the current segment.
+   * Return the quaternion at the given @p time. The quaternion is interpolated
+   * linearly (Slerp) between the start and the end of the segment.
    */
-  Quaternion get_current_quaternion() const;
+  Quaternion get_quaternion(double const time) const;
 
   /**
    * Return the scan path's list of segments
@@ -138,7 +139,8 @@ private:
   void update_current_segment_info(double time, bool save_segment,
                                    unsigned int &current_segment,
                                    dealii::Point<3> &segment_start_point,
-                                   double &segment_start_time) const;
+                                   double &segment_start_time,
+                                   Quaternion &segment_start_rotation) const;
 
   /**
    * Flag is true if we have reached the end of _scan_path_file.

tests/CMakeLists.txt

@@ -116,6 +116,7 @@ adamantine_COPY_DATA_FILE(raytracing_experimental_data_0_0.csv)
 adamantine_COPY_DATA_FILE(raytracing_non_AA_cells-0-0.csv)
 adamantine_COPY_DATA_FILE(scan_path.txt)
 adamantine_COPY_DATA_FILE(scan_path_5_axis.txt)
+adamantine_COPY_DATA_FILE(scan_path_arcs.txt)
 adamantine_COPY_DATA_FILE(scan_path_diagonal.txt)
 adamantine_COPY_DATA_FILE(scan_path_event_series.inp)
 adamantine_COPY_DATA_FILE(scan_path_L.txt)

tests/data/scan_path_arcs.txt

@@ -0,0 +1,4 @@
+0.1, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0
+1.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0
+2.0, 0.5, 0.0, 0.0, 0.70710678118, 0.0, 0.70710678118, 0.0, 1.0
+3.0, 1.0, 0.0, 0.0, 0.70710678118, 0.70710678118, 0.0, 0.0, 1.0

tests/test_scan_path.cc

@@ -1,4 +1,4 @@
-/* SPDX-FileCopyrightText: Copyright (c) 2016 - 2021, the adamantine authors.
+/* SPDX-FileCopyrightText: Copyright (c) 2016 - 2026, the adamantine authors.
  * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
  */
 
@@ -103,4 +103,53 @@ BOOST_AUTO_TEST_CASE(scan_path_location, *utf::tolerance(1e-10))
   BOOST_TEST(power == 0.0);
 }
 
+BOOST_AUTO_TEST_CASE(five_axis, *utf::tolerance(1e-9))
+{
+  boost::optional<boost::property_tree::ptree const &> units_optional_database;
+  ScanPath scan_path("scan_path_arcs.txt", "event_series",
+                     units_optional_database);
+
+  double time = 0.5;
+  auto p = scan_path.value(time);
+  auto rotated_p = scan_path.rotate(time, p);
+  BOOST_TEST(rotated_p[0] == 0.2222222222);
+  BOOST_TEST(rotated_p[1] == 0.0);
+  BOOST_TEST(rotated_p[2] == 0.0);
+
+  time = 1.0;
+  p = scan_path.value(time);
+  rotated_p = scan_path.rotate(time, p);
+  BOOST_TEST(rotated_p[0] == 0.5);
+  BOOST_TEST(rotated_p[1] == 0.0);
+  BOOST_TEST(rotated_p[2] == 0.0);
+
+  time = 1.5;
+  p = scan_path.value(time);
+  rotated_p = scan_path.rotate(time, p);
+  BOOST_TEST(rotated_p[0] == 0.35355339059327373);
+  BOOST_TEST(rotated_p[1] == 0.0);
+  BOOST_TEST(rotated_p[2] == -0.35355339059327373);
+
+  time = 2.0;
+  p = scan_path.value(time);
+  rotated_p = scan_path.rotate(time, p);
+  BOOST_TEST(rotated_p[0] == 0.0);
+  BOOST_TEST(rotated_p[1] == 0.0);
+  BOOST_TEST(rotated_p[2] == -0.5);
+
+  time = 2.5;
+  p = scan_path.value(time);
+  rotated_p = scan_path.rotate(time, p);
+  BOOST_TEST(rotated_p[0] == 0.5);
+  BOOST_TEST(rotated_p[1] == 0.25);
+  BOOST_TEST(rotated_p[2] == -0.5);
+
+  time = 3.0;
+  p = scan_path.value(time);
+  rotated_p = scan_path.rotate(time, p);
+  BOOST_TEST(rotated_p[0] == 1.0);
+  BOOST_TEST(rotated_p[1] == 0.0);
+  BOOST_TEST(rotated_p[2] == 0.0);
+}
+
 } // namespace adamantine