interface_area.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_INTERFACE_AREA_HPP
#define PECLET_FLOW_VOF_INTERFACE_AREA_HPP
 
#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
 
#include "vof/curvature.hpp"     // core: hfColumnHeight, curvFrame, PvFit, paraboloidKappa, ...
#include "vof/phase_change.hpp"  // plicArea
#include "vof/plic.hpp"
 
namespace peclet::flow::vof {
 
enum InterfaceAreaMode : int {
  kAreaPlic = 0,       
  kAreaMetric = 1,     
  kAreaNormal = 2,     
  kAreaFootprint = 3,  
  // ---- WO-P3d: a JOINED surface (marching tetrahedra on the cell-centre lattice). No per-cell
  // construction converges (see `vof/marching_cubes.hpp`); these four do not build a per-cell
  // piece at all, they cut ONE watertight sheet and book each triangle to a cell.
  kAreaMcColour = 4,       
  kAreaMcColourSplit = 5,  
  kAreaMcPlic = 6,         
  kAreaMcPlicSplit = 7,    
};
 
KOKKOS_INLINE_FUNCTION double hfAreaElement(const double h[9]) {
  const double hx = 0.5 * (h[2 + 3 * 1] - h[0 + 3 * 1]);
  const double hy = 0.5 * (h[1 + 3 * 2] - h[1 + 3 * 0]);
  return Kokkos::sqrt(1.0 + hx * hx + hy * hy);
}
 
KOKKOS_INLINE_FUNCTION double hfAreaElement(const double h[9], const VofMetric& g, int d) {
  const int d1 = (d + 1) % 3, d2 = (d + 2) % 3;
  const double hx = 0.5 * (h[2 + 3 * 1] - h[0 + 3 * 1]) * (g.h[d] / g.h[d1]);
  const double hy = 0.5 * (h[1 + 3 * 2] - h[1 + 3 * 0]) * (g.h[d] / g.h[d2]);
  return Kokkos::sqrt(1.0 + hx * hx + hy * hy);
}
 
KOKKOS_INLINE_FUNCTION double hfFootprintArea(const double h[9]) {
  const double hx = 0.5 * (h[2 + 3 * 1] - h[0 + 3 * 1]);
  const double hy = 0.5 * (h[1 + 3 * 2] - h[1 + 3 * 0]);
  const double h0 = h[1 + 3 * 1];
  const double off = 0.5 * (hx + hy) - h0;
  const double f = plicVolume(hx, hy, 0.0, 0.5 + off) - plicVolume(hx, hy, 0.0, -0.5 + off);
  return f * Kokkos::sqrt(1.0 + hx * hx + hy * hy);
}
 
KOKKOS_INLINE_FUNCTION double hfFootprintArea(const double h[9], const VofMetric& g, int d) {
  const int d1 = (d + 1) % 3, d2 = (d + 2) % 3;
  const double hx = 0.5 * (h[2 + 3 * 1] - h[0 + 3 * 1]);
  const double hy = 0.5 * (h[1 + 3 * 2] - h[1 + 3 * 0]);
  const double h0 = h[1 + 3 * 1];
  const double off = 0.5 * (hx + hy) - h0;
  const double f = plicVolume(hx, hy, 0.0, 0.5 + off) - plicVolume(hx, hy, 0.0, -0.5 + off);
  return f * (g.h[d1] * g.h[d2]) * hfAreaElement(h, g, d);
}
 
KOKKOS_INLINE_FUNCTION void hfSurfaceNormal(const double h[9], int d, double m[3]) {
  const double hx = 0.5 * (h[2 + 3 * 1] - h[0 + 3 * 1]);
  const double hy = 0.5 * (h[1 + 3 * 2] - h[1 + 3 * 0]);
  const int d1 = (d + 1) % 3, d2 = (d + 2) % 3;
  const double inv = 1.0 / Kokkos::sqrt(1.0 + hx * hx + hy * hy);
  m[d1] = -hx * inv;
  m[d2] = -hy * inv;
  m[d] = inv;
}
 
KOKKOS_INLINE_FUNCTION void hfSurfaceNormal(const double h[9], int d, const VofMetric& g,
                                            double m[3]) {
  const int d1 = (d + 1) % 3, d2 = (d + 2) % 3;
  const double hx = 0.5 * (h[2 + 3 * 1] - h[0 + 3 * 1]) * (g.h[d] / g.h[d1]);
  const double hy = 0.5 * (h[1 + 3 * 2] - h[1 + 3 * 0]) * (g.h[d] / g.h[d2]);
  const double inv = 1.0 / Kokkos::sqrt(1.0 + hx * hx + hy * hy);
  m[d1] = -hx * inv;
  m[d2] = -hy * inv;
  m[d] = inv;
}
 
KOKKOS_INLINE_FUNCTION void pvSurfaceNormal(const double a[6], const double t1[3],
                                            const double t2[3], const double nn[3], double m[3]) {
  const double inv = 1.0 / Kokkos::sqrt(1.0 + a[1] * a[1] + a[2] * a[2]);
  for (int i = 0; i < 3; ++i)
    m[i] = (-a[1] * t1[i] - a[2] * t2[i] + nn[i]) * inv;
}
 
KOKKOS_INLINE_FUNCTION double interfaceAreaFromNormal(int mode, double mx, double my, double mz,
                                                      double alpha, double c, const double ns[3],
                                                      const VofMetric& g) {
  const double m[3] = {mx, my, mz};
  const double aPlic = plicAreaMetric(mx, my, mz, alpha, g);
  if (mode == kAreaPlic)
    return aPlic;
  if (mode == kAreaNormal) {
    double mi[3];
    vofIndexNormal(ns, g, mi);  // the accurate PHYSICAL normal, pushed to index space
    const double al = plicAlpha(mi[0], mi[1], mi[2], c);
    return plicAreaMetric(mi[0], mi[1], mi[2], al, g);
  }
  // kAreaMetric: keep the PLIC footprint, replace the slope.
  int d = 0;
  double best = Kokkos::fabs(ns[0]);
  if (Kokkos::fabs(ns[1]) > best) {
    best = Kokkos::fabs(ns[1]);
    d = 1;
  }
  if (Kokkos::fabs(ns[2]) > best) {
    best = Kokkos::fabs(ns[2]);
    d = 2;
  }
  if (!(best > 0.0))
    return aPlic;
  double np[3];
  if (!(vofPhysNormal(m, g, np) > 0.0))
    return aPlic;
  return aPlic * Kokkos::fabs(np[d]) / best;
}
 
KOKKOS_INLINE_FUNCTION double interfaceAreaFromNormal(int mode, double mx, double my, double mz,
                                                      double alpha, double c, const double ns[3]) {
  return interfaceAreaFromNormal(mode, mx, my, mz, alpha, c, ns, VofMetric{});
}
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_INTERFACE_AREA_HPP