WIP: framework for gaussian heat source

source/CMakeLists.txt

@@ -6,6 +6,7 @@ set(Adamantine_HEADERS
     ${CMAKE_CURRENT_SOURCE_DIR}/DataAssimilator.hh
     ${CMAKE_CURRENT_SOURCE_DIR}/ElectronBeamHeatSource.hh
     ${CMAKE_CURRENT_SOURCE_DIR}/ExperimentalData.hh
+    ${CMAKE_CURRENT_SOURCE_DIR}/GaussianHeatSource.hh
     ${CMAKE_CURRENT_SOURCE_DIR}/Geometry.hh
     ${CMAKE_CURRENT_SOURCE_DIR}/GoldakHeatSource.hh
     ${CMAKE_CURRENT_SOURCE_DIR}/HeatSource.hh
@@ -42,6 +43,7 @@ set(Adamantine_SOURCES
     ${CMAKE_CURRENT_SOURCE_DIR}/CubeHeatSource.cc
     ${CMAKE_CURRENT_SOURCE_DIR}/DataAssimilator.cc
     ${CMAKE_CURRENT_SOURCE_DIR}/ElectronBeamHeatSource.cc
+    ${CMAKE_CURRENT_SOURCE_DIR}/GaussianHeatSource.cc
     ${CMAKE_CURRENT_SOURCE_DIR}/Geometry.cc
     ${CMAKE_CURRENT_SOURCE_DIR}/GoldakHeatSource.cc
     ${CMAKE_CURRENT_SOURCE_DIR}/MaterialPropertyInstDev.cc

source/GaussianHeatSource.cc

@@ -0,0 +1,211 @@
+/* SPDX-FileCopyrightText: Copyright (c)  2026, the adamantine authors.
+ * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+ */
+
+#include <GaussianHeatSource.hh>
+#include <instantiation.hh>
+#include <types.hh>
+
+#include <deal.II/base/utilities.h>
+
+namespace adamantine
+{
+
+template <int dim>
+GaussianHeatSource<dim>::GaussianHeatSource(
+    boost::property_tree::ptree const &beam_database,
+    boost::optional<boost::property_tree::ptree const &> const
+        &units_optional_database)
+    : HeatSource<dim>(beam_database, units_optional_database),
+      _five_axis(this->_scan_path.is_five_axis())
+{
+  _A = beam_database.get<double>("A");
+  _B = beam_database.get<double>("B");
+}
+
+template <int dim>
+void GaussianHeatSource<dim>::update_time(double time)
+{
+  double const segment_power_modifier =
+      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();
+
+  for (unsigned int d = 0; d < 3; ++d)
+  {
+    _beam_center(d) = path(d);
+  }
+
+  _depth = this->_beam.depth;
+  _radius_squared = this->_beam.radius_squared;
+
+  double const aspect_ratio = std::max(1.0, _depth[0] / this->_beam.radius);
+  double const n =
+      std::min(std::max(0.0, _A * std::log2(aspect_ratio) + _B), 9.0);
+  _k = std::pow(2.0, n);
+
+  double const V0 =
+      0.5 * dealii::numbers::PI * _radius_squared[0] * _depth[0] *
+      std::tgamma(1.0 / _k) / (_k * std::pow(3.0, 1.0 / _k));
+
+  double const effective_power =
+      this->_beam.absorption_efficiency * this->_beam.max_power *
+      segment_power_modifier;
+
+  _alpha = effective_power / V0;
+}
+
+template <int dim>
+double GaussianHeatSource<dim>::value(dealii::Point<dim> const &point) const
+{
+  dealii::Point<dim> rotated_point;
+  if constexpr (dim == 2)
+  {
+    rotated_point = point;
+  }
+  else
+  {
+    rotated_point = _five_axis ? _quaternion.rotate(point) : point;
+  }
+
+  double const z = rotated_point[axis<dim>::z] - _beam_center[axis<dim>::z][0];
+
+  if (z > 0. || z < -_depth[0])
+  {
+    return 0.;
+  }
+
+  double xpy_squared = dealii::Utilities::fixed_power<2>(
+      rotated_point[axis<dim>::x] - _beam_center[axis<dim>::x][0]);
+  if constexpr (dim == 3)
+  {
+    xpy_squared += dealii::Utilities::fixed_power<2>(
+        rotated_point[axis<dim>::y] - _beam_center[axis<dim>::y][0]);
+  }
+
+  if (xpy_squared > 5. * _radius_squared[0])
+  {
+    return 0.;
+  }
+
+  double const radial_component =
+      std::exp(-2.0 * xpy_squared / _radius_squared[0]);
+
+  double const depth_component =
+      std::exp(-3.0 * std::pow(std::abs(z / _depth[0]), _k));
+
+  return _alpha[0] * radial_component * depth_component;
+}
+
+template <int dim>
+dealii::VectorizedArray<double> GaussianHeatSource<dim>::value(
+    dealii::Point<dim, dealii::VectorizedArray<double>> const &points) const
+{
+  dealii::Point<dim, dealii::VectorizedArray<double>> rotated_points;
+  if constexpr (dim == 2)
+  {
+    rotated_points = points;
+  }
+  else
+  {
+    rotated_points = _five_axis ? _quaternion.rotate(points) : points;
+  }
+
+  auto const z = rotated_points[axis<dim>::z] - _beam_center[axis<dim>::z];
+
+  dealii::VectorizedArray<double> depth_mask;
+  for (unsigned int i = 0; i < depth_mask.size(); ++i)
+  {
+    depth_mask[i] = (z[i] > 0. || z[i] < -_depth[i]) ? 0. : 1.;
+  }
+  if (depth_mask.sum() == 0)
+  {
+    return depth_mask;
+  }
+
+  auto xpy_squared = dealii::Utilities::fixed_power<2>(
+      rotated_points[axis<dim>::x] - _beam_center[axis<dim>::x]);
+  if constexpr (dim == 3)
+  {
+    auto const y_squared = dealii::Utilities::fixed_power<2>(
+        rotated_points[axis<dim>::y] - _beam_center[axis<dim>::y]);
+    xpy_squared += y_squared;
+  }
+
+  auto const xpy_area = 5. * _radius_squared - xpy_squared;
+  dealii::VectorizedArray<double> xpy_mask;
+  for (unsigned int i = 0; i < xpy_mask.size(); ++i)
+  {
+    xpy_mask[i] = xpy_area[i] < 0. ? 0. : 1.;
+  }
+  if (xpy_mask.sum() == 0)
+  {
+    return xpy_mask;
+  }
+
+  dealii::VectorizedArray<double> shape = _k;
+
+  auto const radial_component =
+      std::exp(-2.0 * xpy_squared / _radius_squared);
+
+  auto const depth_component =
+      std::exp(-3.0 * std::pow(std::abs(z / _depth), shape));
+
+  return depth_mask * _alpha * radial_component * depth_component;
+}
+
+template <int dim>
+dealii::BoundingBox<dim>
+GaussianHeatSource<dim>::get_bounding_box(double const time,
+                                          double const scaling_factor) const
+{
+  dealii::Point<3> const &beam_center = this->_scan_path.value(time, false);
+  if constexpr (dim == 2)
+  {
+    return {{{beam_center[axis<dim>::x] - scaling_factor * this->_beam.radius,
+              beam_center[axis<dim>::z] - scaling_factor * this->_beam.depth},
+             {beam_center[axis<dim>::x] + scaling_factor * this->_beam.radius,
+              beam_center[axis<dim>::z]}}};
+  }
+  else
+  {
+    dealii::Point<3> max_corner(
+        beam_center[axis<dim>::x] + scaling_factor * this->_beam.radius,
+        beam_center[axis<dim>::y] + scaling_factor * this->_beam.radius,
+        beam_center[axis<dim>::z]);
+    dealii::Point<3> min_corner(
+        beam_center[axis<dim>::x] - scaling_factor * this->_beam.radius,
+        beam_center[axis<dim>::y] - scaling_factor * this->_beam.radius,
+        beam_center[axis<dim>::z] - scaling_factor * this->_beam.depth);
+
+    if (this->_scan_path.is_five_axis())
+    {
+      dealii::Point<3> rotated_max_corner =
+          this->_scan_path.rotate(time, max_corner);
+      dealii::Point<3> rotated_min_corner =
+          this->_scan_path.rotate(time, min_corner);
+
+      dealii::Point<3> new_max_corner;
+      dealii::Point<3> new_min_corner;
+      for (int d = 0; d < 3; ++d)
+      {
+        new_max_corner[d] =
+            std::max(rotated_max_corner[d], rotated_min_corner[d]);
+        new_min_corner[d] =
+            std::min(rotated_max_corner[d], rotated_min_corner[d]);
+      }
+
+      return {{new_min_corner, new_max_corner}};
+    }
+    else
+    {
+      return {{min_corner, max_corner}};
+    }
+  }
+}
+
+} // namespace adamantine
+
+INSTANTIATE_DIM(GaussianHeatSource)

source/GaussianHeatSource.hh

@@ -0,0 +1,87 @@
+/* SPDX-FileCopyrightText: Copyright (c) 2026, the adamantine authors.
+ * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+ */
+
+#ifndef GAUSSIAN_HEAT_SOURCE_HH
+#define GAUSSIAN_HEAT_SOURCE_HH
+
+#include <HeatSource.hh>
+#include <Quaternion.hh>
+
+#include <limits>
+
+namespace adamantine
+{
+/**
+ * A derived class from HeatSource for the Gaussian model of a laser heat source.
+ */
+template <int dim>
+class GaussianHeatSource final : public HeatSource<dim>
+{
+public:
+  /**
+   * Constructor.
+   * \param[in] beam_database requires the following entries:
+   *   - <B>absorption_efficiency</B>: double in \f$[0,1]\f$
+   *   - <B>depth</B>: double in \f$[0,\infty)\f$
+   *   - <B>diameter</B>: double in \f$[0,\infty)\f$
+   *   - <B>max_power</B>: double in \f$[0, \infty)\f$
+   *   - <B>scan_path_file</B>: name of the file that contains the scan path
+   *     segments
+   *   - <B>scan_path_file_format</B>: format of the scan path file
+   *   - <B>A</B>: double, empirical coefficient for the model
+   *   - <B>B</B>: double, empirical coefficient for the model
+   * \param[in] units_optional_database may contain the following entries:
+   *   - <B>heat_source.dimension</B>
+   *   - <B>heat_source.power</B>
+   */
+  GaussianHeatSource(
+      boost::property_tree::ptree const &beam_database,
+      boost::optional<boost::property_tree::ptree const &> const
+          &units_optional_database);
+
+  /**
+   * Set the time variable. This updates the beam position and recalculates
+   * the shape factor (_k) and normalization constant (_alpha).
+   */
+  void update_time(double time) final;
+
+  /**
+   * Returns the value of a Gaussian heat source at a specified point.
+   */
+  double value(dealii::Point<dim> const &point) const final;
+
+  /**
+   * Same function as above but it uses vectorized data.
+   */
+  dealii::VectorizedArray<double>
+  value(dealii::Point<dim, dealii::VectorizedArray<double>> const &points)
+      const final;
+
+  /**
+   * Returns the bounding box of the heat source at a given time.
+   */
+  dealii::BoundingBox<dim>
+  get_bounding_box(double time, double const scaling_factor) const final;
+
+private:
+  double _A;
+  double _B;
+
+  bool const _five_axis;
+  Quaternion _quaternion;
+
+  dealii::Point<3, dealii::VectorizedArray<double>> _beam_center;
+  dealii::VectorizedArray<double> _depth =
+      std::numeric_limits<double>::signaling_NaN();
+  dealii::VectorizedArray<double> _radius_squared =
+      std::numeric_limits<double>::signaling_NaN();
+
+  dealii::VectorizedArray<double> _alpha =
+      std::numeric_limits<double>::signaling_NaN();
+
+  double _k = std::numeric_limits<double>::signaling_NaN();
+};
+} // namespace adamantine
+
+#endif

source/ThermalPhysics.templates.hh

@@ -5,6 +5,7 @@
 #ifndef THERMAL_PHYSICS_TEMPLATES_HH
 #define THERMAL_PHYSICS_TEMPLATES_HH
 
+#include <GaussianHeatSource.hh>
 #include <CubeHeatSource.hh>
 #include <ElectronBeamHeatSource.hh>
 #include <GoldakHeatSource.hh>
@@ -402,6 +403,11 @@ ThermalPhysics<dim, n_materials, p_order, fe_degree, MaterialStates,
       _heat_sources[i] = std::make_shared<CubeHeatSource<dim>>(
           beam_database, units_optional_database);
     }
+    else if (type == "gaussian")
+    {
+      _heat_sources[i] = std::make_shared<GaussianHeatSource<dim>>(
+          beam_database, units_optional_database);
+    }
     else
     {
       ASSERT_THROW(false, "Beam type '" +

source/ensemble_management.cc

@@ -2,6 +2,7 @@
  * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
  */
 
+#include <GaussianHeatSource.hh>
 #include <CubeHeatSource.hh>
 #include <ElectronBeamHeatSource.hh>
 #include <GoldakHeatSource.hh>
@@ -158,6 +159,16 @@ std::vector<std::shared_ptr<HeatSource<dim>>> get_bounding_heat_sources(
       bounding_heat_sources[i] = std::make_shared<CubeHeatSource<dim>>(
           beam_database, units_optional_database);
     }
+    else if (type == "gaussian")
+    {
+      boost::property_tree::ptree modified_beam_database = beam_database;
+      modified_beam_database.put("diameter", diameter_max[bounding_source]);
+      modified_beam_database.put("depth", depth_max[bounding_source]);
+      ++bounding_source;
+
+      bounding_heat_sources[i] = std::make_shared<GaussianHeatSource<dim>>(
+          modified_beam_database, units_optional_database);
+    }
   }
 
   return bounding_heat_sources;

source/validate_input_database.cc

@@ -369,10 +369,11 @@ void validate_input_database(boost::property_tree::ptree &database)
         "sources.beam_" + std::to_string(beam_index) + ".type");
     ASSERT_THROW(boost::iequals(beam_type, "goldak") ||
                      boost::iequals(beam_type, "electron_beam") ||
-                     boost::iequals(beam_type, "cube"),
+                     boost::iequals(beam_type, "cube") ||
+                     boost::iequals(beam_type, "gaussian"),
                  "Beam type, '" + beam_type +
                      "', is not recognized. Valid options are: 'goldak', "
-                     "'electron_beam', and 'cube'.");
+                     "'electron_beam', 'cube', and 'gaussian'.");
     ASSERT_THROW(database.get_child("sources")
                          .get_child("beam_" + std::to_string(beam_index))
                          .count("scan_path_file") != 0,