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