dune-grid-glue 2.8.0
Loading...
Searching...
No Matches
projection_impl.hh
Go to the documentation of this file.
1#include <dune/common/fmatrix.hh>
2
3#include <cmath>
4
5namespace Dune {
6namespace GridGlue {
7
8namespace ProjectionImplementation {
9
20template<typename Coordinate, typename Field>
21inline Coordinate
22corner(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
40inline std::pair<unsigned, unsigned>
41edgeToCorners(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
66template<typename Coordinate, typename Corners>
67inline typename Corners::value_type
68interpolate(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
87template<typename Coordinate, typename Normals>
88inline typename Normals::value_type
89interpolate_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
107template<typename Coordinate, typename Field>
108inline bool
109inside(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
126template<typename Coordinate>
127Projection<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
135template<typename Coordinate>
136void
138::epsilon(const Field epsilon)
139{
140 m_epsilon = epsilon;
141}
142
143template<typename Coordinate>
144template<typename Corners, typename Normals>
145void
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
235template<typename Coordinate>
236template<typename Corners, typename Normals>
237void
238Projection<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
305template<typename Coordinate>
306template<typename Corners, typename Normals>
307void
308Projection<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
433template<typename Coordinate>
434template<typename Corners, typename Normals>
435bool 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
464template<typename Coordinate>
465template<typename Corners, typename Normals>
466void Projection<Coordinate>
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
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
std::pair< unsigned, unsigned > edgeToCorners(unsigned edge)
Definition: projection_impl.hh:41
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:59