signed_distance_solver.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) 2025-2026 The Lethe Authors
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 
#pragma once
#include <deal.II/base/bounding_box.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/table_handler.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/utilities.h>
 
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
 
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/mapping_q.h>
 
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
 
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/q_collection.h>
 
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/lapack_full_matrix.h>
#include <deal.II/lac/vector.h>
 
#include <deal.II/matrix_free/fe_point_evaluation.h>
 
#include <deal.II/non_matching/fe_immersed_values.h>
#include <deal.II/non_matching/fe_values.h>
#include <deal.II/non_matching/mesh_classifier.h>
#include <deal.II/non_matching/quadrature_generator.h>
 
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
 
#include <algorithm>
#include <cmath>
#include <fstream>
#include <iostream>
#include <map>
#include <memory>
#include <optional>
#include <set>
#include <string>
#include <vector>
 
namespace MeltPoolDG::LevelSet
{
  enum class Verbosity
  {
    quiet,
    verbose
  };
 
  template <typename T>
  [[nodiscard]] constexpr int
  sgn(const T val) noexcept
  {
    return (static_cast<T>(0) < val) - (val < static_cast<T>(0));
  }
 
  inline void
  announce_string(const dealii::ConditionalOStream &pcout,
                  const std::string                &expression,
                  const char                        delimiter = '-')
  {
    pcout << std::string(expression.size() + 1, delimiter) << std::endl;
    pcout << expression << std::endl;
    pcout << std::string(expression.size() + 1, delimiter) << std::endl;
  }
 
  template <int dim>
  inline double
  compute_cell_diameter(const double cell_measure, const unsigned int fe_degree)
  {
    if constexpr (dim == 2)
      return std::sqrt(4. * cell_measure / dealii::numbers::PI) / fe_degree;
    else if constexpr (dim == 3)
      return std::cbrt(6. * cell_measure / dealii::numbers::PI) / fe_degree;
    else
      AssertThrow(false, dealii::ExcMessage("Only dim=2 and dim=3 are supported."));
    return 0.0;
  }
 
  template <int dim>
  void
  vertices_cell_mapping(
    const dealii::DoFHandler<dim> &dof_handler,
    std::map<unsigned int, std::set<typename dealii::DoFHandler<dim>::active_cell_iterator>>
      &vertices_cell_map)
  {
    vertices_cell_map.clear();
    for (const auto &cell : dof_handler.active_cell_iterators())
      if (cell->is_locally_owned() || cell->is_ghost())
        for (unsigned int i = 0; i < dealii::GeometryInfo<dim>::vertices_per_cell; ++i)
          vertices_cell_map[cell->vertex_index(i)].insert(cell);
  }
 
  template <int dim>
  std::vector<typename dealii::DoFHandler<dim>::active_cell_iterator>
  find_cells_around_cell(
    std::map<unsigned int, std::set<typename dealii::DoFHandler<dim>::active_cell_iterator>>
                                                                 &vertices_cell_map,
    const typename dealii::DoFHandler<dim>::active_cell_iterator &cell)
  {
    std::set<typename dealii::DoFHandler<dim>::active_cell_iterator> neighbors_cells;
    for (unsigned int i = 0; i < dealii::GeometryInfo<dim>::vertices_per_cell; ++i)
      {
        const unsigned int v_index = cell->vertex_index(i);
        neighbors_cells.insert(vertices_cell_map[v_index].begin(),
                               vertices_cell_map[v_index].end());
      }
 
    return {neighbors_cells.begin(), neighbors_cells.end()};
  }
 
  template <int dim>
  double
  find_point_triangle_distance(const std::vector<dealii::Point<dim>> &triangle,
                               const dealii::Point<dim>              &point)
  {
    double distance = 0.0;
 
    const dealii::Point<dim> &point_0 = triangle[0];
    const dealii::Point<dim> &point_1 = triangle[1];
 
    if constexpr (dim == 3)
      {
        const dealii::Point<dim> &point_2 = triangle[2];
 
        const dealii::Tensor<1, dim> e_0 = point_1 - point_0;
        const dealii::Tensor<1, dim> e_1 = point_2 - point_0;
 
        const double a   = e_0.norm_square();
        const double b   = scalar_product(e_0, e_1);
        const double c   = e_1.norm_square();
        const double det = a * c - b * b;
 
        const dealii::Tensor<1, dim> vector_to_plane = point_0 - point;
 
        const double d = scalar_product(e_0, vector_to_plane);
        const double e = scalar_product(e_1, vector_to_plane);
 
        double s = b * e - c * d;
        double t = b * d - a * e;
 
        if (s + t <= det)
          {
            if (s < 0)
              {
                if (t < 0)
                  {
                    if (d < 0)
                      {
                        t = 0;
                        if (-d >= a)
                          s = 1;
                        else
                          s = -d / a;
                      }
                    else
                      {
                        s = 0;
                        if (e >= 0)
                          t = 0;
                        else if (-e >= c)
                          t = 1;
                        else
                          t = e / c;
                      }
                  }
                else
                  {
                    s = 0;
                    if (e >= 0)
                      t = 0;
                    else if (-e >= c)
                      t = 1;
                    else
                      t = -e / c;
                  }
              }
            else if (t < 0)
              {
                t = 0;
                if (d >= 0)
                  s = 0;
                else if (-d >= a)
                  s = 1;
                else
                  s = -d / a;
              }
            else
              {
                const double inv_det = 1. / det;
                s *= inv_det;
                t *= inv_det;
              }
          }
        else
          {
            if (s < 0)
              {
                const double tmp0 = b + d;
                const double tmp1 = c + e;
                if (tmp1 > tmp0)
                  {
                    const double numer = tmp1 - tmp0;
                    const double denom = a - 2 * b + c;
                    if (numer >= denom)
                      s = 1;
                    else
                      s = numer / denom;
                    t = 1 - s;
                  }
                else
                  {
                    s = 0;
                    if (tmp1 <= 0)
                      t = 1;
                    else if (e >= 0)
                      t = 0;
                    else
                      t = -e / c;
                  }
              }
            else if (t < 0)
              {
                const double tmp0 = b + e;
                const double tmp1 = a + d;
                if (tmp1 > tmp0)
                  {
                    const double numer = tmp1 - tmp0;
                    const double denom = a - 2 * b + c;
                    if (numer >= denom)
                      t = 1;
                    else
                      t = numer / denom;
                    s = 1 - t;
                  }
                else
                  {
                    t = 0;
                    if (tmp1 <= 0)
                      s = 1;
                    else if (d >= 0)
                      s = 0;
                    else
                      s = -d / a;
                  }
              }
            else
              {
                const double numer = (c + e) - (b + d);
                if (numer <= 0)
                  s = 0;
                else
                  {
                    const double denom = a - 2 * b + c;
                    if (numer >= denom)
                      s = 1;
                    else
                      s = numer / denom;
                  }
                t = 1 - s;
              }
          }
 
        const dealii::Point<dim> pt_in_triangle = point_0 + s * e_0 + t * e_1;
        distance                                = pt_in_triangle.distance(point);
      }
 
    if constexpr (dim == 2)
      {
        const dealii::Tensor<1, dim> d     = point_1 - point_0;
        const double                 t_bar = d * (point - point_0) / (d.norm() * d.norm() + 1e-12);
 
        if (t_bar <= 0.0)
          distance = (point - point_0).norm();
        else if (t_bar >= 1.0)
          distance = (point - point_1).norm();
        else
          distance = (point - (point_0 + t_bar * d)).norm();
      }
 
    return distance;
  }
 
  template <int dim,
            typename VectorType = dealii::Vector<double>,
            typename FEType     = dealii::FE_Q<dim>>
  class CellWiseFunction : public dealii::Function<dim>
  {
  public:
    CellWiseFunction(const unsigned int p_fe_degree)
      : fe(p_fe_degree)
      , n_cell_wise_dofs(fe.dofs_per_cell)
    {}
 
    void
    set_active_cell(const VectorType &in_local_dof_values)
    {
      cell_dof_values = in_local_dof_values;
    }
 
    double
    value(const dealii::Point<dim> &point, const unsigned int component = 0) const override
    {
      double value = 0.;
      for (unsigned int i = 0; i < n_cell_wise_dofs; ++i)
        value += cell_dof_values[i] * fe.shape_value_component(i, point, component);
      return value;
    }
 
    dealii::Tensor<1, dim>
    gradient(const dealii::Point<dim> &point, const unsigned int component = 0) const override
    {
      dealii::Tensor<1, dim> gradient;
      for (unsigned int i = 0; i < n_cell_wise_dofs; ++i)
        gradient += cell_dof_values[i] * fe.shape_grad_component(i, point, component);
      return gradient;
    }
 
    dealii::SymmetricTensor<2, dim>
    hessian(const dealii::Point<dim> &point, const unsigned int component = 0) const override
    {
      dealii::Tensor<2, dim> hessian;
      for (unsigned int i = 0; i < n_cell_wise_dofs; ++i)
        hessian += cell_dof_values[i] * fe.shape_grad_grad_component(i, point, component);
      return symmetrize(hessian);
    }
 
  private:
    FEType fe;
    unsigned int n_cell_wise_dofs;
    VectorType cell_dof_values;
  };
 
  template <int dim>
  double
  compute_cell_wise_volume(dealii::FEPointEvaluation<1, dim> &fe_point_evaluation,
                           const typename dealii::DoFHandler<dim>::active_cell_iterator &cell,
                           dealii::Vector<double> cell_dof_level_set_values,
                           const double           corr,
                           const unsigned int     n_quad_points);
 
  template <int dim, typename VectorType>
  std::pair<double, double>
  integrate_volume_and_surface(const dealii::DoFHandler<dim>    &dof_handler,
                               const dealii::FiniteElement<dim> &fe,
                               const VectorType                 &level_set_vector_relevant_copy);
 
  template <int dim, typename VectorType>
  std::pair<double, double>
  integrate_volume_and_surface(const dealii::DoFHandler<dim>    &dof_handler,
                               const dealii::FiniteElement<dim> &fe,
                               const VectorType                 &level_set_vector,
                               const double                      iso_level);
 
  template <int dim, typename VectorType>
  void
  reconstruct_interface(
    const dealii::Mapping<dim>       &mapping,
    const dealii::DoFHandler<dim>    &dof_handler,
    const dealii::FiniteElement<dim> &fe,
    const VectorType                 &level_set_vector,
    const double                      iso_level,
    const int                         mca_n_subdivisions,
    std::map<dealii::types::global_cell_index, std::vector<dealii::Point<dim>>>
      &interface_reconstruction_vertices,
    std::map<dealii::types::global_cell_index,
             std::vector<dealii::CellData<dim == 1 ? 1 : dim - 1>>> &interface_reconstruction_cells,
    std::set<dealii::types::global_dof_index>                       &intersected_dofs);
 
  template <int dim>
  class InterfaceReconstructionDataOut : public dealii::DataOutInterface<0, dim>
  {
  public:
    void
    build_patches(const std::map<dealii::types::global_cell_index, std::vector<dealii::Point<dim>>>
                    &interface_reconstruction_vertices)
    {
      for (const auto &cell : interface_reconstruction_vertices)
        for (const dealii::Point<dim> &vertex : cell.second)
          {
            dealii::DataOutBase::Patch<0, dim> patch;
            patch.vertices[0] = vertex;
            patches.push_back(patch);
          }
    }
 
  private:
    const std::vector<dealii::DataOutBase::Patch<0, dim>> &
    get_patches() const override
    {
      return patches;
    }
 
    std::vector<std::string>
    get_dataset_names() const override
    {
      return dataset_names;
    }
 
    std::vector<dealii::DataOutBase::Patch<0, dim>> patches;
    std::vector<std::string> dataset_names;
  };
 
  template <int dim, typename VectorType = dealii::LinearAlgebra::distributed::Vector<double>>
  class SignedDistanceSolver
  {
  public:
    SignedDistanceSolver(
      std::shared_ptr<dealii::parallel::DistributedTriangulationBase<dim>> background_triangulation,
      const dealii::FiniteElement<dim>                                    &background_fe,
      const double                                                         p_max_distance,
      const double                                                         p_iso_level,
      const double                                                         p_scaling = 1.0,
      const Verbosity p_verbosity = Verbosity::quiet);
 
    SignedDistanceSolver(const dealii::DoFHandler<dim> &background_dof_handler,
                         const double                   p_max_distance,
                         const double                   p_iso_level,
                         const double                   p_scaling   = 1.0,
                         const Verbosity                p_verbosity = Verbosity::quiet);
 
    void
    setup_dofs();
 
    void
    set_level_set_from_background_mesh(const dealii::DoFHandler<dim> &background_dof_handler,
                                       const VectorType              &background_level_set_vector);
 
    void
    solve();
 
    VectorType &
    get_signed_distance();
 
    const VectorType &
    get_signed_distance() const;
 
    const dealii::DoFHandler<dim> &
    get_dof_handler() const
    {
      return *dof_handler;
    }
 
    void
    output_interface_reconstruction(const std::string &filename) const;
 
  private:
    void
    zero_out_ghost_values() const;
 
    void
    update_ghost_values() const;
 
    void
    exchange_distance();
 
    void
    initialize_distance();
 
    void
    compute_first_neighbors_distance();
 
    void
    compute_second_neighbors_distance();
 
    void
    compute_signed_distance_from_distance();
 
    void
    compute_cell_wise_volume_correction();
 
    void
    conserve_global_volume();
 
    inline void
    set_face_opposite_dofs_map();
 
    inline void
    set_face_dofs_map();
 
    inline void
    get_face_opposite_dofs(unsigned int               local_face_id,
                           std::vector<unsigned int> &local_opposite_dofs) const;
 
    inline void
    get_face_local_dofs(unsigned int local_face_id, std::vector<unsigned int> &local_dofs) const;
 
    inline void
    get_face_transformation_jacobian(
      const dealii::DerivativeForm<1, dim, dim> &cell_transformation_jac,
      const unsigned int                         local_face_id,
      dealii::LAPACKFullMatrix<double>          &face_transformation_jac) const;
 
    inline dealii::Point<dim>
    transform_ref_face_point_to_ref_cell(const dealii::Point<dim - 1> &x_ref_face,
                                         const unsigned int            local_face_id) const;
 
    inline void
    compute_residual(const dealii::Tensor<1, dim>           &x_n_to_x_I_real,
                     const dealii::Tensor<1, dim>           &distance_gradient,
                     const dealii::LAPACKFullMatrix<double> &transformation_jac,
                     dealii::Tensor<1, dim - 1>             &residual_ref) const;
 
    inline dealii::Tensor<1, dim>
    transform_ref_face_correction_to_ref_cell(const dealii::Vector<double> &correction_ref_face,
                                              const unsigned int            local_face_id) const;
 
    inline void
    compute_analytical_jacobian(const dealii::Tensor<1, dim>           &x_n_to_x_I_real_p1,
                                const dealii::LAPACKFullMatrix<double> &transformation_jacobian,
                                const std::vector<double>              &face_local_dof_values,
                                dealii::LAPACKFullMatrix<double>       &jacobian_matrix);
 
    inline double
    compute_distance(const dealii::Tensor<1, dim> &x_n_to_x_I_real, const double distance) const
    {
      return distance + x_n_to_x_I_real.norm();
    }
    bool is_external_dof_handler = false;
    std::unique_ptr<dealii::FiniteElement<dim>> fe;
    std::unique_ptr<dealii::DoFHandler<dim>> owned_dof_handler;
 
    const dealii::DoFHandler<dim> *dof_handler = nullptr;
    std::unique_ptr<dealii::MappingQ<dim>> mapping;
 
    const double max_distance;
    const double iso_level;
    const double scaling;
    Verbosity verbosity;
 
    dealii::ConditionalOStream pcout;
    dealii::IndexSet locally_owned_dofs;
    dealii::IndexSet locally_relevant_dofs;
    dealii::IndexSet locally_active_dofs;
 
    VectorType level_set;
    VectorType signed_distance;
    VectorType signed_distance_with_ghost;
    VectorType distance;
    VectorType distance_with_ghost;
    VectorType volume_correction;
 
    std::map<dealii::types::global_cell_index, std::vector<dealii::Point<dim>>>
      interface_reconstruction_vertices;
    std::map<dealii::types::global_cell_index,
             std::vector<dealii::CellData<dim == 1 ? 1 : dim - 1>>>
      interface_reconstruction_cells;
    std::set<dealii::types::global_dof_index> intersected_dofs;
 
    std::map<unsigned int, std::vector<unsigned int>> face_opposite_dofs_map;
    std::map<unsigned int, std::vector<unsigned int>> face_dofs_map;
 
    std::unique_ptr<dealii::AffineConstraints<double>> constraints;
  };
} // namespace MeltPoolDG::LevelSet