distance_functions.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>
 
#include <algorithm>
#include <array>
#include <cmath>
 
 
namespace MeltPoolDG::Functions::Distance
{
  template <int dim, typename number>
  class Plane : public dealii::Function<dim, number>
  {
  public:
    Plane(const dealii::Point<dim, number>     &center_of_wall,
          const dealii::Tensor<1, dim, number> &normal,
          const std::array<number, dim - 1>    &half_extents,
          const dealii::Tensor<1, dim, number> &first_tangent_direction =
            dealii::Tensor<1, dim, number>())
      : center_of_wall(center_of_wall)
      , normal(normal)
      , half_extents(half_extents)
    {
      AssertThrow(std::abs(normal.norm() - 1.) < tolerance,
                  dealii::ExcMessage("The normal of the wall plane must be a unit vector."));
 
      // Building  tangent basis
      if constexpr (dim == 2)
        {
          tangents[0][0] = -normal[1];
          tangents[0][1] = normal[0];
        }
      else if constexpr (dim == 3)
        {
          const number first_tangent_norm = first_tangent_direction.norm();
          AssertThrow(first_tangent_norm > tolerance,
                      dealii::ExcMessage("The first tangent direction must be nonzero in 3D."));
 
          tangents[0] = first_tangent_direction / first_tangent_norm;
          tangents[1] = dealii::cross_product_3d(normal, tangents[0]);
          AssertThrow(std::abs(normal * tangents[0]) < tolerance,
                      dealii::ExcMessage(
                        "The first tangent direction must be perpendicular to the normal."));
        }
    }
 
    dealii::Tensor<1, dim, number>
    gradient(const dealii::Point<dim, number> &x, const unsigned int component = 0) const override
    {
      Assert(component == 0,
             dealii::ExcMessage("The distance function is scalar; component must be 0."));
 
      const dealii::Tensor<1, dim, number> r      = offset_from_wall(x);
      const number                         r_norm = r.norm();
 
      // At the wall surface, r vanishes and r/r_norm is undefined. Fall back to the wall's unit
      // normal as the gradient direction.
      if (r_norm < tolerance)
        return normal;
 
      return r / r_norm;
    }
 
    number
    value(const dealii::Point<dim, number> &x, const unsigned int component = 0) const override
    {
      Assert(component == 0,
             dealii::ExcMessage("The distance function is scalar; component must be 0."));
 
      return offset_from_wall(x).norm();
    }
 
  private:
    dealii::Tensor<1, dim, number>
    offset_from_wall(const dealii::Point<dim, number> &x) const
    {
      const dealii::Tensor<1, dim, number> position_relative_to_center = x - center_of_wall;
      const number                         normal_coord = position_relative_to_center * normal;
 
      dealii::Tensor<1, dim, number> offset = normal_coord * normal;
      for (unsigned int i = 0; i < dim - 1; ++i)
        {
          const number tangential_coord = position_relative_to_center * tangents[i];
          const number clamped_tangential_coord =
            std::clamp(tangential_coord, -half_extents[i], half_extents[i]);
          offset += (tangential_coord - clamped_tangential_coord) * tangents[i];
        }
 
      return offset;
    }
 
    dealii::Point<dim, number> center_of_wall;
 
    dealii::Tensor<1, dim, number> normal;
 
    std::array<dealii::Tensor<1, dim, number>, dim - 1> tangents;
 
    std::array<number, dim - 1> half_extents;
 
    const number tolerance = 1e-12;
  };
 
  template <int dim, typename number>
  class TruncatedCone : public dealii::Function<dim, number>
  {
  public:
    TruncatedCone(const dealii::Point<dim, number>     &axis_center,
                  const dealii::Tensor<1, dim, number> &axis_direction,
                  const number                          radius_bottom,
                  const number                          radius_top,
                  const number                          half_length)
      : axis_center(axis_center)
      , radius_bottom(radius_bottom)
      , radius_top(radius_top)
      , half_length(half_length)
    {
      AssertThrow(axis_direction.norm() > tolerance,
                  dealii::ExcMessage("The axis direction of the truncated cone must not be zero."));
      AssertThrow(half_length > 0.,
                  dealii::ExcMessage("The half length of the truncated cone must be positive."));
      AssertThrow(radius_bottom >= 0. && radius_top >= 0.,
                  dealii::ExcMessage("The radii of the truncated cone must not be negative."));
      AssertThrow(radius_bottom > 0. || radius_top > 0.,
                  dealii::ExcMessage(
                    "At least one radius of the truncated cone must be greater than zero."));
 
      axis = axis_direction / axis_direction.norm(); // unit vector
    }
 
    dealii::Tensor<1, dim, number>
    gradient(const dealii::Point<dim, number> &x, const unsigned int component = 0) const override
    {
      Assert(component == 0,
             dealii::ExcMessage("The distance function is scalar; component must be 0."));
 
      const dealii::Tensor<1, dim, number> d =
        x - axis_center; // Vector from the axis center to the point
      const dealii::Tensor<1, dim, number> d_radial =
        d - (d * axis) * axis; // Vector from the axis (perpendicular) to the point.
      const number radial_coord = d_radial.norm();
 
      dealii::Tensor<1, dim, number> radial_direction;
      if (radial_coord >= tolerance)
        radial_direction = d_radial / radial_coord;
      else
        {
          // The Cartesian unit vector least aligned with the axis, made perpendicular to it is
          // used, so the gradient stays well-defined.
          unsigned int least_aligned = 0;
          for (unsigned int i = 1; i < dim; ++i)
            if (std::abs(axis[i]) < std::abs(axis[least_aligned]))
              least_aligned = i;
 
          dealii::Tensor<1, dim, number> unit_vector;
          unit_vector[least_aligned] = 1.;
          radial_direction           = unit_vector - (unit_vector * axis) * axis;
          radial_direction /= radial_direction.norm();
        }
 
      const auto [axial_offset, radial_offset] = meridian_offset_from_wall(x);
      const number distance                    = std::hypot(axial_offset, radial_offset);
 
      // On the surface, fall back to the inward surface normal.
      if (distance < tolerance)
        {
          const number axial_span   = 2. * half_length;
          const number radial_span  = radius_top - radius_bottom;
          const number slant_length = std::hypot(axial_span, radial_span);
 
          return (radial_span * axis - axial_span * radial_direction) / slant_length;
        }
 
      return (radial_offset * radial_direction + axial_offset * axis) / distance;
    }
 
    number
    value(const dealii::Point<dim, number> &x, const unsigned int component = 0) const override
    {
      Assert(component == 0,
             dealii::ExcMessage("The distance function is scalar; component must be 0."));
 
      const auto [axial_offset, radial_offset] = meridian_offset_from_wall(x);
      return std::hypot(axial_offset, radial_offset);
    }
 
  private:
    std::array<number, 2>
    meridian_offset_from_wall(const dealii::Point<dim, number> &x) const
    {
      const dealii::Tensor<1, dim, number> d =
        x - axis_center;                   // Vector from the axis center to the point
      const number axial_coord = d * axis; // Coordinate of the point along the axis.
      const number radial_coord =
        (d - axial_coord * axis).norm(); // Coordinate of the point from the axis.
 
      const number axial_span  = 2. * half_length;
      const number radial_span = radius_top - radius_bottom;
 
      // Project (axial_coord, radial_coord) onto the meridian segment, clamped to its end points.
      const number segment_fraction =
        std::clamp(((axial_coord + half_length) * axial_span +
                    (radial_coord - radius_bottom) * radial_span) /
                     (axial_span * axial_span + radial_span * radial_span),
                   number(0.),
                   number(1.));
 
      const number axial_coord_closest  = -half_length + segment_fraction * axial_span;
      const number radial_coord_closest = radius_bottom + segment_fraction * radial_span;
 
      return {{axial_coord - axial_coord_closest, radial_coord - radial_coord_closest}};
    }
 
    const dealii::Point<dim, number> axis_center;
 
    dealii::Tensor<1, dim, number> axis;
 
    const number radius_bottom;
 
    const number radius_top;
 
    const number half_length;
 
    const number tolerance = 1e-12;
  };
} // namespace MeltPoolDG::Functions::Distance