dune-grid-glue  2.8.0
projection_impl.hh
Go to the documentation of this file.
1 #include <dune/common/fmatrix.hh>
2 
3 #include <cmath>
4 
5 namespace Dune {
6 namespace GridGlue {
7 
8 namespace ProjectionImplementation {
9 
20 template<typename Coordinate, typename Field>
21 inline Coordinate
22 corner(unsigned c)
23 {
24  Coordinate x(Field(0));
25  if (c == 0)
26  return x;
27  x[c-1] = Field(1);
28  return x;
29 }
30 
40 inline std::pair<unsigned, unsigned>
41 edgeToCorners(unsigned edge)
42 {
43  switch(edge) {
44  case 0: return {0, 1};
45  case 1: return {0, 2};
46  case 2: return {1, 2};
47  }
48  DUNE_THROW(Dune::Exception, "Unexpected edge number.");
49 }
50 
66 template<typename Coordinate, typename Corners>
67 inline typename Corners::value_type
68 interpolate(const Coordinate& x, const Corners& corners)
69 {
70  auto y = corners[0];
71  for (unsigned i = 0; i < corners.size() - 1; ++i)
72  y.axpy(x[i], corners[i+1] - corners[0]);
73  return y;
74 }
75 
87 template<typename Coordinate, typename Normals>
88 inline typename Normals::value_type
89 interpolate_unit_normals(const Coordinate& x, const Normals& normals)
90 {
91  auto n = interpolate(x, normals);
92  n /= n.two_norm();
93  return n;
94 }
95 
107 template<typename Coordinate, typename Field>
108 inline bool
109 inside(const Coordinate& x, const Field& epsilon)
110 {
111  const unsigned dim = Coordinate::dimension;
112  Field sum(0);
113  for (unsigned i = 0; i < dim-1; ++i) {
114  if (x[i] < -epsilon)
115  return false;
116  sum += x[i];
117  }
118  /* If any xᵢ is NaN, sum will be NaN and this comparison false! */
119  if (sum <= Field(1) + epsilon)
120  return true;
121  return false;
122 }
123 
124 } /* namespace ProjectionImplementation */
125 
126 template<typename Coordinate>
128 ::Projection(const Field overlap, const Field max_normal_product)
129  : m_overlap(overlap)
130  , m_max_normal_product(max_normal_product)
131 {
132  /* Nothing. */
133 }
134 
135 template<typename Coordinate>
136 void
138 ::epsilon(const Field epsilon)
139 {
140  m_epsilon = epsilon;
141 }
142 
143 template<typename Coordinate>
144 template<typename Corners, typename Normals>
145 void
147 ::doProjection(const std::tuple<Corners&, Corners&>& corners, const std::tuple<Normals&, Normals&>& normals)
148 {
149  /* Try to obtain Φ(xᵢ) for each corner xᵢ of the preimage triangle.
150  * This means solving a linear system of equations
151  * Φ(xᵢ) = (1-α-β) y₀ + α y₁ + β y₂ = xᵢ + δ nᵢ
152  * or α (y₁ - y₀) + β (y₂ - y₀) - δ nᵢ = xᵢ - y₀
153  * to obtain the barycentric coordinates (α, β) of Φ(xᵢ) in the image
154  * triangle and the distance δ.
155  *
156  * In the matrix m corresponding to the system, only the third column and the
157  * right-hand side depend on i. The first two columns can be assembled before
158  * and reused.
159  */
160  using namespace ProjectionImplementation;
161  using std::get;
162  typedef Dune::FieldMatrix<Field, dim, dim> Matrix;
163  Matrix m;
164 
165  const auto& origin = get<0>(corners);
166  const auto& origin_normals = get<0>(normals);
167  const auto& target = get<1>(corners);
168  const auto& target_normals = get<1>(normals);
169  auto& images = get<0>(m_images);
170  auto& success = get<0>(m_success);
171 
172  /* directionsᵢ = (yᵢ - y₀) / ||yᵢ - y₀||
173  * These are the first to columns of the system matrix; the rescaling is done
174  * to ensure all columns have a comparable norm (the last has the normal with norm 1.
175  */
176  std::array<Coordinate, dim-1> directions;
177  std::array<Field, dim-1> scales;
178  /* estimator for the diameter of the target face */
179  Field scaleSum(0);
180  for (unsigned i = 0; i < dim-1; ++i) {
181  directions[i] = target[i+1] - target[0];
182  scales[i] = directions[i].infinity_norm();
183  directions[i] /= scales[i];
184  scaleSum += scales[i];
185  }
186 
187  for (unsigned i = 0; i < dim-1; ++i) {
188  for (unsigned j = 0; j < dim; ++j) {
189  m[j][i] = directions[i][j];
190  }
191  }
192 
193  m_projection_valid = true;
194  success.reset();
195 
196  /* Now project xᵢ for each i */
197  for (unsigned i = 0; i < origin.size(); ++i) {
198  for (unsigned j = 0; j < dim; ++j)
199  m[j][dim-1] = origin_normals[i][j];
200 
201  const Coordinate rhs = origin[i] - target[0];
202 
203  try {
204  /* y = (α, β, δ) */
205  auto& y = images[i];
206  m.solve(y, rhs);
207  for (unsigned j = 0; j < dim-1; ++j)
208  y[j] /= scales[j];
209  /* Solving gave us -δ as the term is "-δ nᵢ". */
210  y[dim-1] *= Field(-1);
211 
212  /* If the forward projection is too far in the wrong direction
213  * then this might result in artificial inverse projections or
214  * edge intersections. To prevent these wrong cases but not
215  * dismiss feasible intersections, the projection is dismissed
216  * if the forward projection is further than two times the
217  * approximate diameter of the image triangle.
218  */
219  if(y[dim-1] < -2*scaleSum) {
220  success.set(i,false);
221  m_projection_valid = false;
222  return;
223  }
224 
225  const bool feasible = projectionFeasible(origin[i], origin_normals[i], y, target, target_normals);
226  success.set(i, feasible);
227  }
228  catch (const Dune::FMatrixError&) {
229  success.set(i, false);
230  m_projection_valid = false;
231  }
232  }
233 }
234 
235 template<typename Coordinate>
236 template<typename Corners, typename Normals>
237 void
238 Projection<Coordinate>
239 ::doInverseProjection(const std::tuple<Corners&, Corners&>& corners, const std::tuple<Normals&, Normals&>& normals)
240 {
241  /* Try to obtain Φ⁻¹(yᵢ) for each corner yᵢ of the image triangle.
242  * Instead of solving the problem directly (which would lead to
243  * non-linear equations), we make use of the forward projection Φ
244  * which projects the preimage triangle on the plane spanned by the
245  * image triangle. The inverse projection is then given by finding
246  * the barycentric coordinates of yᵢ with respect to the triangle
247  * with the corners Φ(xᵢ). This way we only have to solve linear
248  * equations.
249  */
250 
251  using namespace ProjectionImplementation;
252  using std::get;
253  typedef Dune::FieldMatrix<Field, dim-1, dim-1> Matrix;
254  typedef Dune::FieldVector<Field, dim-1> Vector;
255 
256  /* The inverse projection can only be computed if the forward projection
257  * managed to project all xᵢ on the plane spanned by the yᵢ
258  */
259  if (!m_projection_valid) {
260  get<1>(m_success).reset();
261  return;
262  }
263 
264  const auto& images = get<0>(m_images);
265  const auto& target_corners = get<1>(corners);
266  auto& preimages = get<1>(m_images);
267  auto& success = get<1>(m_success);
268 
269  std::array<Coordinate, dim> v;
270  for (unsigned i = 0; i < dim-1; ++i) {
271  v[i] = interpolate(images[i+1], target_corners);
272  v[i] -= interpolate(images[0], target_corners);
273  }
274 
275  Matrix m;
276  for (unsigned i = 0; i < dim-1; ++i) {
277  for (unsigned j = 0; j < dim-1; ++j) {
278  m[i][j] = v[i]*v[j];
279  }
280  }
281 
282  for (unsigned i = 0; i < dim; ++i) {
283  /* Convert yᵢ to barycentric coordinates with respect to Φ(xⱼ) */
284  v[dim-1] = target_corners[i];
285  v[dim-1] -= interpolate(images[0], target_corners);
286 
287  Vector rhs, z;
288  for (unsigned j = 0; j < dim-1; ++j)
289  rhs[j] = v[dim-1]*v[j];
290  m.solve(z, rhs);
291 
292  for (unsigned j = 0; j < dim-1; ++j)
293  preimages[i][j] = z[j];
294 
295  /* Calculate distance along normal direction */
296  const auto x = interpolate(z, get<0>(corners));
297  preimages[i][dim-1] = (x - target_corners[i]) * get<1>(normals)[i];
298 
299  /* Check y_i lies inside the Φ(xⱼ) */
300  const bool feasible = projectionFeasible(target_corners[i], get<1>(normals)[i], preimages[i], get<0>(corners), get<0>(normals));
301  success.set(i, feasible);
302  }
303 }
304 
305 template<typename Coordinate>
306 template<typename Corners, typename Normals>
307 void
308 Projection<Coordinate>
309 ::doEdgeIntersection(const std::tuple<Corners&, Corners&>& corners, const std::tuple<Normals&, Normals&>& normals)
310 {
311  using namespace ProjectionImplementation;
312  using std::get;
313 
314  m_number_of_edge_intersections = 0;
315 
316  /* There are no edge intersections for 2d, only for 3d */
317  if (dim != 3)
318  return;
319 
320  /* There are no edge intersections
321  * - when the projection is invalid,
322  * - when the projected triangle lies fully in the target triangle,
323  * - or when the target triangle lies fully in the projected triangle.
324  */
325  if (!m_projection_valid || get<0>(m_success).all() || get<1>(m_success).all()) {
326  return;
327  }
328 
329  const auto& images = get<0>(m_images);
330  const auto& ys = get<1>(corners);
331 
332  /* Intersect line through Φ(xᵢ), Φ(xⱼ) with line through yₖ, yₗ:
333  We want α, β ∈ ℝ such that
334  Φ(xᵢ) + α (Φ(xⱼ) - Φ(xᵢ)) = yₖ + β (yₗ - yₖ)
335  or
336  α (Φ(xⱼ)-Φ(xᵢ)) + β (yₗ-yₖ) = yₖ-Φ(xᵢ)
337  To get a 2×2 system of equations, multiply with yₘ-y₀ for
338  m ∈ {1,̣̣2} which are linear indep. (and so the system is
339  equivalent to the original 3×2 system)
340  */
341  for (unsigned edgex = 0; edgex < dim; ++edgex) {
342  unsigned i, j;
343  std::tie(i, j) = edgeToCorners(edgex);
344 
345  /* Both sides of edgex lie in the target triangle means no edge intersection */
346  if (get<0>(m_success)[i] && get<0>(m_success)[j])
347  continue;
348 
349  const auto pxi = interpolate(images[i], ys);
350  const auto pxj = interpolate(images[j], ys);
351  const auto pxjpxi = pxj - pxi;
352 
353  typedef Dune::FieldMatrix<Field, dim-1, dim-1> Matrix;
354  typedef Dune::FieldVector<Field, dim-1> Vector;
355 
356  for (unsigned edgey = 0; edgey < dim; ++edgey) {
357  unsigned k, l;
358  std::tie(k, l) = edgeToCorners(edgey);
359 
360  /* Both sides of edgey lie in the projected triangle means no edge intersection */
361  if (get<1>(m_success)[k] && get<1>(m_success)[l])
362  continue;
363 
364  const auto ykyl = ys[k] - ys[l];
365  const auto ykpxi = ys[k] - pxi;
366 
367  /* If edges are parallel then the intersection is already computed by vertex projections. */
368  bool parallel = true;
369  for (unsigned h=0; h<3; h++)
370  parallel &= std::abs(ykyl[(h+1)%3]*pxjpxi[(h+2)%3] - ykyl[(h+2)%3]*pxjpxi[(h+1)%3])<1e-14;
371  if (parallel)
372  continue;
373 
374  Matrix mat;
375  Vector rhs, z;
376 
377  for (unsigned m = 0; m < dim-1; ++m) {
378  const auto ym1y0 = ys[m+1] - ys[0];
379  mat[m][0] = pxjpxi * ym1y0;
380  mat[m][1] = ykyl * ym1y0;
381  rhs[m] = ykpxi * ym1y0;
382  }
383 
384  try {
385  using std::isfinite;
386 
387  mat.solve(z, rhs);
388 
389  /* If solving the system gives a NaN, the edges are probably parallel. */
390  if (!isfinite(z[0]) || !isfinite(z[1]))
391  continue;
392 
393  /* Filter out corner (pre)images. We only want "real" edge-edge intersections here. */
394  if (z[0] < m_epsilon || z[0] > Field(1) - m_epsilon
395  || z[1] < m_epsilon || z[1] > Field(1) - m_epsilon)
396  continue;
397 
398  Coordinate local_x = corner<Coordinate, Field>(i);
399  local_x.axpy(z[0], corner<Coordinate, Field>(j) - corner<Coordinate, Field>(i));
400  Coordinate local_y = corner<Coordinate, Field>(k);
401  local_y.axpy(z[1], corner<Coordinate, Field>(l) - corner<Coordinate, Field>(k));
402 
403  /* Make sure the intersection is in the triangle. */
404  if (!inside(local_x, m_epsilon) || !inside(local_y, m_epsilon))
405  continue;
406 
407  /* Make sure the intersection respects overlap. */
408  auto xy = interpolate(local_x, get<0>(corners));
409  xy -= interpolate(local_y, get<1>(corners));
410  const auto nx = interpolate_unit_normals(local_x, get<0>(normals));
411  const auto ny = interpolate_unit_normals(local_y, get<1>(normals));
412  local_x[dim-1] = -(xy*nx);
413  local_y[dim-1] = xy*ny;
414 
415  if (local_x[dim-1] < -m_overlap-m_epsilon || local_y[dim-1] < -m_overlap-m_epsilon)
416  continue;
417 
418  /* Normals should be opposing. */
419  if (nx*ny > m_max_normal_product + m_epsilon)
420  continue;
421 
422  /* Intersection is feasible. Store it. */
423  auto& intersection = m_edge_intersections[m_number_of_edge_intersections++];
424  intersection = { {{edgex, edgey}}, {{local_x, local_y}} };
425  }
426  catch(const Dune::FMatrixError&) {
427  /* Edges might be parallel, ignore and continue with next edge */
428  }
429  }
430  }
431 }
432 
433 template<typename Coordinate>
434 template<typename Corners, typename Normals>
435 bool Projection<Coordinate>
436 ::projectionFeasible(const Coordinate& x, const Coordinate& nx, const Coordinate& px, const Corners& corners, const Normals& normals) const
437 {
438  using namespace ProjectionImplementation;
439 
440  /* Image must be within simplex. */
441  if (!inside(px, m_epsilon))
442  return false;
443 
444  /* Distance along normal must not be smaller than -overlap. */
445  if (px[dim-1] < -m_overlap-m_epsilon)
446  return false;
447 
448  /* Distance along normal at image must not be smaller than -overlap. */
449  auto xmy = x;
450  xmy -= interpolate(px, corners);
451  const auto n = interpolate_unit_normals(px, normals);
452  const auto d = xmy * n;
453  if (d < -m_overlap-m_epsilon)
454  return false;
455 
456  /* Normals at x and Φ(x) are opposing. */
457  if (nx * n > m_max_normal_product + m_epsilon)
458  return false;
459 
460  /* Okay, projection is feasible. */
461  return true;
462 }
463 
464 template<typename Coordinate>
465 template<typename Corners, typename Normals>
467 ::project(const std::tuple<Corners&, Corners&>& corners, const std::tuple<Normals&, Normals&>& normals)
468 {
469  doProjection(corners, normals);
470  doInverseProjection(corners, normals);
471  doEdgeIntersection(corners, normals);
472 }
473 
474 } /* namespace GridGlue */
475 } /* namespace Dune */
Definition: gridglue.hh:35
std::pair< unsigned, unsigned > edgeToCorners(unsigned edge)
Definition: projection_impl.hh:41
Corners::value_type interpolate(const Coordinate &x, const Corners &corners)
Definition: projection_impl.hh:68
bool inside(const Coordinate &x, const Field &epsilon)
Definition: projection_impl.hh:109
Coordinate corner(unsigned c)
Definition: projection_impl.hh:22
Normals::value_type interpolate_unit_normals(const Coordinate &x, const Normals &normals)
Definition: projection_impl.hh:89
Projection of a line (triangle) on another line (triangle).
Definition: projection.hh:19
Coordinate::field_type Field
Scalar type.
Definition: projection.hh:54
Projection(const Field overlap=Field(0), const Field max_normal_product=Field(-0.1))
Definition: projection_impl.hh:128
void epsilon(const Field epsilon)
Set epsilon used for floating-point comparisons.
Definition: projection_impl.hh:138
void project(const std::tuple< Corners &, Corners & > &corners, const std::tuple< Normals &, Normals & > &normals)
Do the actual projection.
Definition: projection_impl.hh:467