implicitsurfaceprojection.hh 2.61 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CURVED_SURFACE_GRID_IMPLICIT_SURFACE_PROJECTION_HH
#define DUNE_CURVED_SURFACE_GRID_IMPLICIT_SURFACE_PROJECTION_HH

#include <cmath>
#include <type_traits>

#include <dune/common/parametertree.hh>

namespace Dune {

//! Closest-point projection to surface given by implicit function
/**
 * Surface S is given by zero-levelset of function F.
 * We assume that F is differentiable in order to evaluate normals and to
 * do an iterative projection.
 **/
template< class F >
class ImplicitSurfaceProjection
{
  using Functor = F;
  using DFunctor = std::decay_t<decltype(derivative(std::declval<F>()))>;

public:
  //! Constructor for a given differentiable function with surface = { x : phi(x) == 0 }
  ImplicitSurfaceProjection (const F& phi, int maxIter = 10)
    : maxIter_(maxIter)
    , phi_(phi)
    , dphi_(derivative(phi_))
  {}

  //! Evaluation of the closest-point projection
  /**
   * Simple iterative algorithm proposed by Demlow and Dziuk in
   * AN ADAPTIVE FINITE ELEMENT METHOD FOR THE LAPLACE-BELTRAMI OPERATOR ON IMPLICITLY DEFINED SURFACES
   **/
  template< class Domain >
  Domain operator() (const Domain& x0) const
  {
    using std::sqrt;
    using T = typename Domain::value_type;

    Domain x = x0;
    T tol = sqrt(std::numeric_limits<T>::epsilon());

    auto phi = phi_(x);
    auto sign = phi > 0 ? 1 : phi < 0 ? -1 : 0;
    auto grad_phi = dphi_(x);
    auto grad_phi_norm2 = grad_phi.two_norm2();

    T err = 1;
    for (int i = 0; i < maxIter_ && err > tol; ++i) {
      auto y = x - grad_phi * (phi/grad_phi_norm2);

      auto dist = sign * (y - x0).two_norm();
      auto grad_phi_y = dphi_(y);

      auto& dx = grad_phi_y;
      dx *= -dist/grad_phi_y.two_norm();
      x = x0 + dx;

      phi = phi_(x);
      grad_phi = dphi_(x);
      grad_phi_norm2 = grad_phi.two_norm2();
      err = sqrt(phi*phi / grad_phi_norm2 + (grad_phi / grad_phi_norm2 - dx / dx.two_norm()).two_norm2());
    }

    return x;
  }

  //! Normal vector given by grad(F)/|grad(F)|
  template< class Domain >
  Domain normal (const Domain& x0) const
  {
    auto grad_phi = dphi_(x0);
    return grad_phi / grad_phi.two_norm();
  }

private:
  int maxIter_;
  Functor phi_;
  DFunctor dphi_;
};

//! Construct a grid function representing a projection to an implicit surface
template< class Grid, class F >
auto implicitSurfaceGridFunction (const F& phi, int maxIter = 10)
{
  return analyticGridFunction<Grid>(ImplicitSurfaceProjection<F>{phi, maxIter});
}

} // end namespace Dune

#endif // DUNE_CURVED_SURFACE_GRID_IMPLICIT_SURFACE_PROJECTION_HH