1#include <dune/common/fmatrix.hh>
8namespace ProjectionImplementation {
20template<
typename Coordinate,
typename Field>
24 Coordinate x(Field(0));
40inline std::pair<unsigned, unsigned>
44 case 0:
return {0, 1};
45 case 1:
return {0, 2};
46 case 2:
return {1, 2};
48 DUNE_THROW(Dune::Exception,
"Unexpected edge number.");
66template<
typename Coordinate,
typename Corners>
67inline typename Corners::value_type
71 for (
unsigned i = 0; i < corners.size() - 1; ++i)
72 y.axpy(x[i], corners[i+1] - corners[0]);
87template<
typename Coordinate,
typename Normals>
88inline typename Normals::value_type
107template<
typename Coordinate,
typename Field>
109inside(
const Coordinate& x,
const Field& epsilon)
111 const unsigned dim = Coordinate::dimension;
113 for (
unsigned i = 0; i < dim-1; ++i) {
119 if (sum <= Field(1) + epsilon)
126template<
typename Coordinate>
127Projection<Coordinate>
128::Projection(
const Field overlap,
const Field max_normal_product)
130 , m_max_normal_product(max_normal_product)
135template<
typename Coordinate>
143template<
typename Coordinate>
144template<
typename Corners,
typename Normals>
147::doProjection(
const std::tuple<Corners&, Corners&>& corners,
const std::tuple<Normals&, Normals&>& normals)
160 using namespace ProjectionImplementation;
162 typedef Dune::FieldMatrix<Field, dim, dim> Matrix;
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);
176 std::array<Coordinate, dim-1> directions;
177 std::array<Field, dim-1> scales;
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];
187 for (
unsigned i = 0; i < dim-1; ++i) {
188 for (
unsigned j = 0; j < dim; ++j) {
189 m[j][i] = directions[i][j];
193 m_projection_valid =
true;
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];
201 const Coordinate rhs = origin[i] - target[0];
207 for (
unsigned j = 0; j < dim-1; ++j)
210 y[dim-1] *= Field(-1);
219 if(y[dim-1] < -2*scaleSum) {
220 success.set(i,
false);
221 m_projection_valid =
false;
225 const bool feasible = projectionFeasible(origin[i], origin_normals[i], y, target, target_normals);
226 success.set(i, feasible);
228 catch (
const Dune::FMatrixError&) {
229 success.set(i,
false);
230 m_projection_valid =
false;
235template<
typename Coordinate>
236template<
typename Corners,
typename Normals>
238Projection<Coordinate>
239::doInverseProjection(
const std::tuple<Corners&, Corners&>& corners,
const std::tuple<Normals&, Normals&>& normals)
251 using namespace ProjectionImplementation;
253 typedef Dune::FieldMatrix<Field, dim-1, dim-1> Matrix;
254 typedef Dune::FieldVector<Field, dim-1> Vector;
259 if (!m_projection_valid) {
260 get<1>(m_success).reset();
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);
269 std::array<Coordinate, dim> v;
270 for (
unsigned i = 0; i < dim-1; ++i) {
276 for (
unsigned i = 0; i < dim-1; ++i) {
277 for (
unsigned j = 0; j < dim-1; ++j) {
282 for (
unsigned i = 0; i < dim; ++i) {
284 v[dim-1] = target_corners[i];
285 v[dim-1] -=
interpolate(images[0], target_corners);
288 for (
unsigned j = 0; j < dim-1; ++j)
289 rhs[j] = v[dim-1]*v[j];
292 for (
unsigned j = 0; j < dim-1; ++j)
293 preimages[i][j] = z[j];
297 preimages[i][dim-1] = (x - target_corners[i]) * get<1>(normals)[i];
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);
305template<
typename Coordinate>
306template<
typename Corners,
typename Normals>
308Projection<Coordinate>
309::doEdgeIntersection(
const std::tuple<Corners&, Corners&>& corners,
const std::tuple<Normals&, Normals&>& normals)
311 using namespace ProjectionImplementation;
314 m_number_of_edge_intersections = 0;
325 if (!m_projection_valid || get<0>(m_success).all() || get<1>(m_success).all()) {
329 const auto& images = get<0>(m_images);
330 const auto& ys = get<1>(corners);
341 for (
unsigned edgex = 0; edgex < dim; ++edgex) {
346 if (get<0>(m_success)[i] && get<0>(m_success)[j])
351 const auto pxjpxi = pxj - pxi;
353 typedef Dune::FieldMatrix<Field, dim-1, dim-1> Matrix;
354 typedef Dune::FieldVector<Field, dim-1> Vector;
356 for (
unsigned edgey = 0; edgey < dim; ++edgey) {
361 if (get<1>(m_success)[k] && get<1>(m_success)[l])
364 const auto ykyl = ys[k] - ys[l];
365 const auto ykpxi = ys[k] - pxi;
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;
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;
390 if (!isfinite(z[0]) || !isfinite(z[1]))
394 if (z[0] < m_epsilon || z[0] > Field(1) - m_epsilon
395 || z[1] < m_epsilon || z[1] > Field(1) - m_epsilon)
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));
404 if (!
inside(local_x, m_epsilon) || !
inside(local_y, m_epsilon))
412 local_x[dim-1] = -(xy*nx);
413 local_y[dim-1] = xy*ny;
415 if (local_x[dim-1] < -m_overlap-m_epsilon || local_y[dim-1] < -m_overlap-m_epsilon)
419 if (nx*ny > m_max_normal_product + m_epsilon)
423 auto& intersection = m_edge_intersections[m_number_of_edge_intersections++];
424 intersection = { {{edgex, edgey}}, {{local_x, local_y}} };
426 catch(
const Dune::FMatrixError&) {
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
438 using namespace ProjectionImplementation;
441 if (!
inside(px, m_epsilon))
445 if (px[dim-1] < -m_overlap-m_epsilon)
452 const auto d = xmy * n;
453 if (d < -m_overlap-m_epsilon)
457 if (nx * n > m_max_normal_product + m_epsilon)
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)
469 doProjection(corners, normals);
470 doInverseProjection(corners, normals);
471 doEdgeIntersection(corners, normals);
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