sdf.h

Go to the documentation of this file.

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2019 Adrian Kummerlaender
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public
 *  License along with this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/
 
#ifndef SDF_H
#define SDF_H
 
#include "core/vector.h"
 
#include <algorithm>
#include <functional>
 
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#ifndef M_PI2
#define M_PI2 1.57079632679489661923
#endif
 
namespace olb {
 
namespace sdf {
 
template <typename T>
T mix(T a, T b, T h) any_platform
{
  return b * (1.0 - h) + a * h;
}
 
template <typename T0, typename T1, typename T2>
decltype(T0 {} * T1 {} * T2 {}) clamp(T0 x, T1 a, T2 b) any_platform
{
  if (x < a) {
    return a;
  }
  else if (x > b) {
    return b;
  }
  else {
    return x;
  }
}
 
template <typename T>
Vector<bool, 3> isSymmetric(std::function<T(const Vector<T, 3>&)> sdf,
                            Vector<T, 3> center) any_platform
{
  Vector<bool, 3> isSymmetric;
 
  isSymmetric[0] = sdf(Vector<T, 3> {1., 1., 1.} + center) ==
                       sdf(Vector<T, 3> {-1., 1., 1.} + center) &&
                   sdf(Vector<T, 3> {1., -1., 1.} + center) ==
                       sdf(Vector<T, 3> {-1., -1., 1.} + center) &&
                   sdf(Vector<T, 3> {1., 1., -1.} + center) ==
                       sdf(Vector<T, 3> {-1., 1., -1.} + center) &&
                   sdf(Vector<T, 3> {1., -1., -1.} + center) ==
                       sdf(Vector<T, 3> {-1., -1., -1.} + center);
  isSymmetric[1] = sdf(Vector<T, 3> {1., 1., 1.} + center) ==
                       sdf(Vector<T, 3> {1., -1., 1.} + center) &&
                   sdf(Vector<T, 3> {-1., 1., 1.} + center) ==
                       sdf(Vector<T, 3> {-1., -1., 1.} + center) &&
                   sdf(Vector<T, 3> {1., 1., -1.} + center) ==
                       sdf(Vector<T, 3> {1., -1., -1.} + center) &&
                   sdf(Vector<T, 3> {-1., 1., -1.} + center) ==
                       sdf(Vector<T, 3> {-1., -1., -1.} + center);
  isSymmetric[2] = sdf(Vector<T, 3> {1., 1., 1.} + center) ==
                       sdf(Vector<T, 3> {1., 1., -1.} + center) &&
                   sdf(Vector<T, 3> {-1., 1., 1.} + center) ==
                       sdf(Vector<T, 3> {-1., 1., -1.} + center) &&
                   sdf(Vector<T, 3> {1., -1., 1.} + center) ==
                       sdf(Vector<T, 3> {1., -1., -1.} + center) &&
                   sdf(Vector<T, 3> {-1., -1., 1.} + center) ==
                       sdf(Vector<T, 3> {-1., -1., -1.} + center);
 
  return isSymmetric;
}
 
template <typename T, unsigned D>
T sphere(Vector<T, D> p, T r) any_platform
{
  return norm(p) - r;
}
 
template <typename T>
T box(Vector<T, 2> p, Vector<T, 2> b) any_platform
{
  Vector<T, 2> q = abs(p) - b;
  return norm(maxv(q, Vector<T, 2>(0.0))) +
         util::min(util::max({q[0], q[1]}), T());
}
 
template <typename T>
T box(Vector<T, 3> p, Vector<T, 3> b) any_platform
{
  Vector<T, 3> q = abs(p) - b;
  return norm(maxv(q, Vector<T, 3>(0.0))) +
         util::min(util::max({q[0], q[1], q[2]}), T());
}
 
template <typename T>
T triangle(Vector<T, 2> p, Vector<T, 2> a, Vector<T, 2> b,
           Vector<T, 2> c) any_platform
{
  Vector<T, 2> e0  = b - a;
  Vector<T, 2> v0  = p - a;
  Vector<T, 2> e1  = c - b;
  Vector<T, 2> v1  = p - b;
  Vector<T, 2> e2  = a - c;
  Vector<T, 2> v2  = p - c;
  Vector<T, 2> pq0 = v0 - e0 * clamp((v0 * e0) / (e0 * e0), 0., 1.);
  Vector<T, 2> pq1 = v1 - e1 * clamp((v1 * e1) / (e1 * e1), 0., 1.);
  Vector<T, 2> pq2 = v2 - e2 * clamp((v2 * e2) / (e2 * e2), 0., 1.);
  T            s   = util::sign(e0[0] * e2[1] - e0[1] * e2[0]);
 
  T dx = util::min(util::min((pq0 * pq0), (pq1 * pq1)), (pq2 * pq2));
  T dy = util::min(util::min(s * (v0[0] * e0[1] - v0[1] * e0[0]),
                             s * (v1[0] * e1[1] - v1[1] * e1[0])),
                   s * (v2[0] * e2[1] - v2[1] * e2[0]));
  return -util::sqrt(dx) * util::sign(dy);
}
 
template <typename T>
T cylinder(Vector<T, 3> p, Vector<T, 3> a, Vector<T, 3> ba, T baba,
           T r) any_platform
{
  const Vector<T, 3> pa   = p - a;
  const T            paba = pa * ba;
  const T            x    = norm(pa * baba - ba * paba) - r * baba;
  const T            y    = util::abs(paba - baba * T {0.5}) - baba * T {0.5};
  const T            x2   = x * x;
  const T            y2   = y * y * baba;
  const T            d    = (util::max(x, y) < T {0})
                                ? -util::min(x2, y2)
                                : (((x > T {0}) ? x2 : T {0}) + ((y > T {0}) ? y2 : T {0}));
  return util::sign(d) * util::sqrt(util::abs(d)) / baba;
}
 
template <typename T>
T cylinder(Vector<T, 3> p, Vector<T, 3> a, Vector<T, 3> b, T r) any_platform
{
  Vector<T, 3> ba   = b - a;
  T            baba = ba * ba;
  return sdf::cylinder(p, a, ba, baba, r);
}
 
template <typename T>
T cone(Vector<T, 3> p, Vector<T, 3> a, Vector<T, 3> ba, T baba, T ra,
       T rb) any_platform
{
  T rba  = rb - ra;
  T papa = (p - a) * (p - a);
  T paba = (p - a) * (ba) / baba;
  // x distance to the axis of the cone
  // note: abs() prevents from negative values, e.g. caused by rounding errors
  T x = util::sqrt(util::abs(papa - paba * paba * baba));
  // cax and cay for the distances to the caps
  T cax = util::max(0.0, x - ((paba < 0.5) ? ra : rb));
  T cay = util::abs(paba - 0.5) - 0.5;
  // cbx and cby for the distance to the side wall
  T k   = rba * rba + baba;
  T f   = clamp((rba * (x - ra) + paba * baba) / k, 0.0, 1.0);
  T cbx = x - ra - f * rba;
  T cby = paba - f;
  T s   = (cbx < 0.0 && cay < 0.0) ? -1.0 : 1.0;
  return s * util::sqrt(util::min(cax * cax + cay * cay * baba,
                                  cbx * cbx + cby * cby * baba));
}
 
template <typename T>
T cone(Vector<T, 3> p, Vector<T, 3> a, Vector<T, 3> b, T ra, T rb) any_platform
{
  Vector<T, 3> ba   = b - a;
  T            baba = ba * ba;
  return sdf::cone(p, a, ba, baba, ra, rb);
}
 
template <typename T>
T ellipsoid(Vector<T, 3> p, Vector<T, 3> r) any_platform
{
  const Vector<T, 3> a(p[0] / r[0], p[1] / r[1], p[2] / r[2]);
  T                  k0 = norm(a);
  const Vector<T, 3> r2(r[0] * r[0], r[1] * r[1], r[2] * r[2]);
  const Vector<T, 3> b(p[0] / r2[0], p[1] / r2[1], p[2] / r2[2]);
  T                  k1 = norm(b);
  return (k0 < 1.0) ? (k0 - 1.0) * util::min(util::min(r[0], r[1]), r[2])
                    : k0 * (k0 - 1.0) / k1;
}
 
template <typename T>
T torus(Vector<T, 3> p, Vector<T, 2> t) any_platform
{
  Vector<T, 2> b {p[0], p[2]};
  Vector<T, 2> q {norm(b) - t[0], p[1]};
  return norm(q) - t[1];
}
 
template <typename T>
T solidAngle(Vector<T, 3> p, Vector<T, 2> c, T r) any_platform
{
  Vector<T, 2> q {norm(Vector<T, 2> {p[0], p[2]}), p[1]};
  T            l = norm(q) - r;
  T            m = norm(q - c * clamp(q * c, 0.0, r));
  return util::max(l, m * util::sign(c[1] * q[0] - c[0] * q[1]));
}
 
template <typename T, unsigned D>
Vector<T, D> translate(Vector<T, D> p, Vector<T, D> origin) any_platform
{
  return p - origin;
}
 
template <typename T>
Vector<T, 3> flip(Vector<T, 3> p) any_platform
{
  return {p[1], p[0], p[2]};
}
 
template <typename T>
T subtraction(T a, T b) any_platform
{
  return util::max(-a, b);
}
 
template <typename T>
T unify(T a, T b) any_platform
{
  return util::min(a, b);
}
 
template <typename T>
T intersection(T a, T b) any_platform
{
  return util::max(a, b);
}
 
template <typename T>
T smooth_union(T a, T b, T k) any_platform
{
  T h = clamp(0.5 + 0.5 * (b - a) / k, 0.0, 1.0);
  return mix(a, b, h) - k * h * (1.0 - h);
}
 
template <typename T>
T smooth_subtraction(T a, T b, T k) any_platform
{
  T h = clamp(0.5 - 0.5 * (b + a) / k, 0.0, 1.0);
  return mix(b, -a, h) + k * h * (1.0 - h);
}
 
template <typename T>
T smooth_intersection(T a, T b, T k) any_platform
{
  T h = clamp(0.5 - 0.5 * (b - a) / k, 0.0, 1.0);
  return mix(b, a, h) + k * h * (1.0 - h);
}
 
template <typename T>
T rounding(T a, T r) any_platform
{
  return a - r;
}
 
template <typename T, bool symmetryCheck = true>
T elongation(std::function<T(const Vector<T, 3>&)> sdf, const Vector<T, 3>& p,
             const Vector<T, 3>& h,
             const Vector<T, 3>& center = (T(0))) any_platform
{
  Vector<T, 3> q       = abs(p - center) - h;
  Vector<T, 3> p_elong = p - center;
  // elongation requires symmetry e.g. for elongation in x --> plane of symmetry = YZ-plane
  Vector<bool, 3> sdfIsSymmetric;
  if constexpr (symmetryCheck) {
    sdfIsSymmetric = isSymmetric(sdf, center);
  }
 
  if (h[0] != 0) {
    p_elong[0] = util::max(q[0], 0.);
    if constexpr (symmetryCheck) {
      if (!sdfIsSymmetric[0]) {
        std::cout << "Warning: symmetry in x is not met" << std::endl;
      }
    }
  }
 
  if (h[1] != 0) {
    p_elong[1] = util::max(q[1], 0.);
    if constexpr (symmetryCheck) {
      if (!sdfIsSymmetric[1]) {
        std::cout << "Warning: symmetry in y is not met" << std::endl;
      }
    }
  }
 
  if (h[2] != 0) {
    p_elong[2] = util::max(q[2], 0.);
    if constexpr (symmetryCheck) {
      if (!sdfIsSymmetric[2]) {
        std::cout << "Warning: symmetry in z is not met" << std::endl;
      }
    }
  }
 
  // Second term needed in case of 3D-Elongation if all parts of Vector q are negative
  return sdf(p_elong + center) + util::min(util::max({q[0], q[1], q[2]}), 0.);
}
 
template <typename T, unsigned D>
T scale(std::function<T(const Vector<T, D>&)> sdf, const Vector<T, D>& p, T s,
        const Vector<T, D>& center = (T(0))) any_platform
{
  return sdf((p - center) / s) * s;
}
 
// TODO: Rename function
template <typename T>
T signedDistanceToPorosity(T signedDist, T eps) any_platform
{
  const T d = signedDist + .5 * eps;
  return util::pow(util::cos(M_PI2 * d / eps), 2);
}
 
template <typename T>
bool evalSolidVolumeFraction(T output[], T signedDist, T eps) any_platform
{
  T const halfEps = .5 * eps;
  if (signedDist <= -halfEps) {
    output[0] = 1.;
    return true;
  }
  else if (signedDist < halfEps) {
    output[0] = signedDistanceToPorosity(signedDist, eps);
    return true;
  }
  output[0] = 0.;
  return false;
}
 
} // namespace sdf
 
} // namespace olb
 
#endif