explicitsurface.hh 11 KB
Newer Older
1
2
#ifndef DUNE_CURVEDGRID_EXPLICIT_SURFACE_PROJECTION_HH
#define DUNE_CURVEDGRID_EXPLICIT_SURFACE_PROJECTION_HH
3
4
5

#include <cmath>

6
7
8
#include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
#include <dune/curvedgrid/utility/kdtree.hh>
#include <dune/curvedgrid/utility/vertexcache.hh>
9

10
#if COLLECTOR_BENCHMARK
11
12
13
  #include <dune/common/timer.hh>
#endif

14
namespace Dune {
15
16
17
18
19
20
21
22

/// \brief Closest-point projection to a surface given by a dune-grid
/**
 * Surface S is given by a triangles-only-grid in a 3d-world. The quality
 * of that grid is crucial to prevent artifacts in the projection. I.e. the
 * grid should be very fine-grained and triangles should not have acute
 * angles. A delaunay-triangulation would be perfect.
 **/
23
template <class T = double>
24
25
26
27
28
class ExplicitSurfaceProjection
{
  using Self = ExplicitSurfaceProjection<T>;
  using Domain = FieldVector<T, 3>;

29
30
  struct SurfaceElement
  {
31
    unsigned v0, v1, v2;
32
33
34
    SurfaceElement (unsigned v0, unsigned v1, unsigned v2)
      : v0(v0), v1(v1), v2(v2)
    {}
35
36
37
  };

  //calculates the distance between a given point and the line segment between 'lin0' and 'lin1'.
38
39
40
41
42
  static T distance_point_line_3d (const Domain& point, const Domain& lin0, const Domain& lin1,
                                                                            Domain &hit_vertex)
  {
    using std::sqrt;

43
44
45
46
47
48
49
50
51
52
53
54
    //the (complete) line through 'lin0' and 'lin1' is given by 'l = lin0 + t*g' with g as follows:
    T gx = lin1[0] - lin0[0];
    T gy = lin1[1] - lin0[1];
    T gz = lin1[2] - lin0[2];

    //the normal-vector from 'point' to l is 'v = lin0 + t*g - point' for some t which is determined
    // by 'v*g = 0'
    T t = (-gx * (lin0[0] - point[0]) - gy * (lin0[1] - point[1]) - gz * (lin0[2] - point[2]))
                                                                      / (gx * gx + gy * gy + gz * gz);

    //if the orthogonal projection of 'point' onto l is not on the given line segment then one of the
    //vertices 'lin0' or 'lin1' is the nearest point to 'point'
55
    if (t < 0) {
56
57
58
59
      //'lin0' is nearest
      hit_vertex[0] = lin0[0];
      hit_vertex[1] = lin0[1];
      hit_vertex[2] = lin0[2];
60
    } else if(t > 1) {
61
62
63
64
      //'lin1' is nearest
      hit_vertex[0] = lin1[0];
      hit_vertex[1] = lin1[1];
      hit_vertex[2] = lin1[2];
65
    } else {
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
      //orthogonal projection is nearest
      hit_vertex[0] = lin0[0] + t * gx;
      hit_vertex[1] = lin0[1] + t * gy;
      hit_vertex[2] = lin0[2] + t * gz;
    }

    //calculate distance
    T dx = point[0] - hit_vertex[0];
    T dy = point[1] - hit_vertex[1];
    T dz = point[2] - hit_vertex[2];

    return sqrt(dx * dx + dy * dy + dz * dz);
  }

  //calculates the distance between a given 'point' and the triangle given by 'tri0', 'tri1', 'tri2'.
81
82
83
84
85
86
  static T distance_point_triangle_3d (const Domain& point,
                    const Domain& tri0, const Domain& tri1, const Domain& tri2, Domain& hit_vertex)
  {
    using std::abs;
    using std::sqrt;
    const T epsilon = sqrt(std::numeric_limits<T>::epsilon());
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109

    //triangle lies in plane 'p = tri0 + u*g + v*h', where g and h are defined as follows:
    T gx = tri1[0] - tri0[0];
    T gy = tri1[1] - tri0[1];
    T gz = tri1[2] - tri0[2];
    T hx = tri2[0] - tri0[0];
    T hy = tri2[1] - tri0[1];
    T hz = tri2[2] - tri0[2];

    //a vector from 'point' to p can be described by 'l = tri0 + u*g + v*h - point' and since we
    //search the distance we need such a vector that is orthogonal to p. Thus 'l*g = 0' and 'l*h = 0'.
    //We get a system of two equations of the form:
    //a1*u + b1*v = c1
    //b1*u + b2*v = c2
    //with a1, b1, b2, c1, c2 as follows:
    T b1 = gx * hx + gy * hy + gz * hz;
    T c2 = -hx * (tri0[0] - point[0]) - hy * (tri0[1] - point[1]) - hz * (tri0[2] - point[2]);
    T b2 = hx * hx + hy * hy + hz * hz;
    T c1 = -gx * (tri0[0] - point[0]) - gy * (tri0[1] - point[1]) - gz * (tri0[2] - point[2]);
    T a1 = gx * gx + gy * gy + gz * gz;

    //solve the system:
    T u, v;
110
    if (abs(b1) < epsilon) {
111
112
      v = c2 / b2;
      u = c1 / a1;
113
    } else {
114
115
116
117
118
119
      v = (c1 * b1 - c2 * a1) / (b1 * b1 - b2 * a1);
      u = (c1 - v * b1) / a1;
    }

    //now test whether the orthogonal projection of 'point' onto p lies inside the triangle
    T distance;
120
    if (u >= 0 && v >= 0 && u + v <= 1) {
121
122
123
124
125
126
127
128
      //yes, inside. The length of the normal-vector is the desired distance
      T fx = hit_vertex[0] = tri0[0] + u * gx + v * hx;
      T fy = hit_vertex[1] = tri0[1] + u * gy + v * hy;
      T fz = hit_vertex[2] = tri0[2] + u * gz + v * hz;
      fx -= point[0];
      fy -= point[1];
      fz -= point[2];
      distance = sqrt(fx * fx + fy * fy + fz * fz);
129
    } else {
130
131
132
133
      //no, not inside. The desired distance is the distance to one of the triangle-edges
      distance = distance_point_line_3d(point, tri0, tri1, hit_vertex);
      Domain current_hit_vertex;
      T current_distance = distance_point_line_3d(point, tri1, tri2, current_hit_vertex);
134
      if (current_distance < distance) {
135
136
137
138
        distance = current_distance;
        hit_vertex = current_hit_vertex;
      }
      current_distance = distance_point_line_3d(point, tri2, tri0, current_hit_vertex);
139
      if (current_distance < distance) {
140
141
142
143
144
145
146
147
        distance = current_distance;
        hit_vertex = current_hit_vertex;
      }
    }
    return distance;
  }

public:
148
  // constructor using a pre-loaded Dune-surface-grid
149
  template <class Grid>
150
  ExplicitSurfaceProjection (const Grid& grid, bool cached = true)
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
    : vertices_(grid.size(Grid::dimension))
    , adjacent_elements_(grid.size(Grid::dimension))
    , minCorner_(std::numeric_limits<T>::max())
    , maxCorner_(std::numeric_limits<T>::min())
    , cached_(cached)
  {
    const auto& indexSet = grid.leafGridView().indexSet();

    //collect vertices for kd-tree
    for (const auto& v : vertices(grid.leafGridView()))
    {
      const auto& vertex = v.geometry().corner(0);
      for (int i = 0; i < Grid::dimensionworld; ++i) {
        if (vertex[i] < minCorner_[i]) minCorner_[i] = vertex[i];
        if (vertex[i] > maxCorner_[i]) maxCorner_[i] = vertex[i];
      }
      vertices_[indexSet.index(v)] = vertex;
    }

    //collect elements and element-adjacency for closest-point-tests
    elements_.reserve(grid.size(0));
    unsigned el_index = 0;
    for (const auto& el : elements(grid.leafGridView()))
    {
      unsigned v0 = indexSet.subIndex(el, 0, Grid::dimension);
      unsigned v1 = indexSet.subIndex(el, 1, Grid::dimension);
      unsigned v2 = indexSet.subIndex(el, 2, Grid::dimension);
      elements_.emplace_back(v0, v1, v2);
      adjacent_elements_[v0].push_back(el_index);
      adjacent_elements_[v1].push_back(el_index);
      adjacent_elements_[v2].push_back(el_index);
      ++el_index;
    }

    kdtree_.emplace(vertices_, minCorner_, maxCorner_);
    kdtree_->update();

188
    if (cached_) {
189
190
191
192
193
      Domain c_expansion = (maxCorner_ - minCorner_) * 0.1;
      cache_.emplace(minCorner_ - c_expansion, maxCorner_ + c_expansion);
    }
  }

194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
  ExplicitSurfaceProjection (const Self& other) = delete;
//     : vertices_(other.vertices_)
//     , adjacent_elements_(other.adjacent_elements_)
//     , elements_(other.elements_)
//     , minCorner_(other.minCorner_)
//     , maxCorner_(other.maxCorner_)
//     , cached_(other.cached_)
//     , cache_(other.cache_)
// #if COLLECTOR_BENCHMARK
//     , inputs_(other.inputs_)
// #endif
//   {
//     kdtree_.emplace(vertices_, minCorner_, maxCorner_);
//     kdtree_->update();
//   }

  ExplicitSurfaceProjection (Self&& other) = delete;
//     : vertices_(std::move(other.vertices_))
//     , adjacent_elements_(std::move(other.adjacent_elements_))
//     , elements_(std::move(other.elements_))
//     , minCorner_(other.minCorner_)
//     , maxCorner_(other.maxCorner_)
//     , cached_(other.cached_)
//     , cache_(std::move(other.cache_))
// #if COLLECTOR_BENCHMARK
//     , inputs_(std::move(other.inputs_))
// #endif
//   {
//     // update pointer
//     kdtree_.emplace(vertices_, minCorner_, maxCorner_);
//     kdtree_->update();
//   }
226

227
228
  Self& operator= (const Self&) = delete;
  Self& operator= (Self&&) = delete;
229

230
  void resetCache ()
231
232
233
234
  {
    if(cached_) cache_->clear();
  }

235
  Domain closestPoint (const Domain& x, std::size_t closest_vertex) const
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
  {
    T min_dist = std::numeric_limits<T>::max();
    Domain hit_vertex, best_hit_vertex;
    for(unsigned el_index : adjacent_elements_[closest_vertex]){
      const SurfaceElement &srf_el = elements_[el_index];
      T current_dist = distance_point_triangle_3d(x,
               vertices_[srf_el.v0], vertices_[srf_el.v1], vertices_[srf_el.v2], hit_vertex);
      if(current_dist < min_dist){
        min_dist = current_dist;
        best_hit_vertex = hit_vertex;
      }
    }
    return best_hit_vertex;
  }

  //evaluation of the closest-point projection
252
  Domain operator() (const Domain& x) const
253
254
255
256
257
258
  {
    #if COLLECTOR_BENCHMARK != 0
      inputs_.push_back(x);
    #endif

    Domain result;
259
    if (!cached_ || !cache_->find(x, result)){
260
261
262
263
264
265
266
267
268
      //find closest surface-vertex
      std::vector<std::size_t> outIndices(1);
      std::vector<T> outDistSq(1);
      kdtree_->query(x, 1, outIndices, outDistSq);

      //check surface-elements adjacent to closest surface-vertex for actual closest-point
      result = closestPoint(x, outIndices[0]);

      //cache result
269
      if (cached_) cache_->insert(x, result);
270
271
272
273
274
    }

    return result;
  }

275
#if COLLECTOR_BENCHMARK
276
277
278
  void collectorBenchmark (bool fresh_cache=true) const
  {
    if (cached_ && fresh_cache) cache_->clear();
279
280
    int uncached_cnt = 0;
    Timer t;
281
    for (auto &input : inputs_){
282
      Domain result;
283
      if (!cached_ || !cache_->find(input, result)) {
284
285
286
287
288
289
290
291
292
        //find closest surface-vertex
        std::vector<std::size_t> outIndices(1);
        std::vector<T> outDistSq(1);
        kdtree_->query(input, 1, outIndices, outDistSq);

        //check surface-elements adjacent to closest surface-vertex for actual closest-point
        result = closestPoint(input, outIndices[0]);

        //cache result
293
        if (cached_) cache_->insert(input, result);
294
295
296
297
298
299
300
        ++uncached_cnt;
      }
      input = result;
    }
    std::cout << "surfGF time: " << t.elapsed() << ", count: "<< inputs_.size()
                                 << ", uncached: " << uncached_cnt << std::endl;
  }
301
#endif
302
303
304
305
306
307
308

private:
  std::vector<Domain> vertices_;
  std::vector<std::vector<unsigned>> adjacent_elements_;
  std::vector<SurfaceElement> elements_;
  Domain minCorner_;
  Domain maxCorner_;
309
  std::optional<KDTree<T>> kdtree_;
310
  bool cached_;
311
  mutable std::optional<VertexCache<Domain, T>> cache_;
312
313
314
315
316
  #if COLLECTOR_BENCHMARK != 0
    mutable std::vector<Domain> inputs_;
  #endif
};

317
318
319
320
321
// deduction guide
template <class Grid>
ExplicitSurfaceProjection(Grid const&, bool = false)
  -> ExplicitSurfaceProjection<typename Grid::ctype>;

322
/// \brief construct a grid function representing an explicit surface parametrization
323
324
template <class Grid>
auto explicitSurfaceProjection (Grid& grid, bool cached = true)
325
{
326
  return ExplicitSurfaceProjection<typename Grid::ctype>{grid, cached};
327
328
329
330
}

} // end namespace Dune

331
#endif // DUNE_CURVEDGRID_EXPLICIT_SURFACE_PROJECTION_HH