1 #include <dune/common/fmatrix.hh> 8 namespace ProjectionImplementation {
20 template<
typename Coordinate,
typename Field>
24 Coordinate x(Field(0));
40 inline 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.");
66 template<
typename Coordinate,
typename Corners>
67 inline typename Corners::value_type
71 for (
unsigned i = 0; i < corners.size() - 1; ++i)
72 y.axpy(x[i], corners[i+1] - corners[0]);
87 template<
typename Coordinate,
typename Normals>
88 inline typename Normals::value_type
107 template<
typename Coordinate,
typename Field>
109 inside(
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)
126 template<
typename Coordinate>
130 , m_max_normal_product(max_normal_product)
135 template<
typename Coordinate>
143 template<
typename Coordinate>
144 template<
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;
235 template<
typename Coordinate>
236 template<
typename Corners,
typename Normals>
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);
305 template<
typename Coordinate>
306 template<
typename Corners,
typename Normals>
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;
370 for (
unsigned m = 0; m < dim-1; ++m) {
371 const auto ym1y0 = ys[m+1] - ys[0];
372 mat[m][0] = pxjpxi * ym1y0;
373 mat[m][1] = ykyl * ym1y0;
374 rhs[m] = ykpxi * ym1y0;
383 if (!isfinite(z[0]) || !isfinite(z[1]))
387 if (z[0] < m_epsilon || z[0] > Field(1) - m_epsilon
388 || z[1] < m_epsilon || z[1] > Field(1) - m_epsilon)
391 Coordinate local_x = corner<Coordinate, Field>(i);
392 local_x.axpy(z[0], corner<Coordinate, Field>(j) - corner<Coordinate, Field>(i));
393 Coordinate local_y = corner<Coordinate, Field>(k);
394 local_y.axpy(z[1], corner<Coordinate, Field>(l) - corner<Coordinate, Field>(k));
397 if (!
inside(local_x, m_epsilon) || !
inside(local_y, m_epsilon))
405 local_x[dim-1] = -(xy*nx);
406 local_y[dim-1] = xy*ny;
408 if (local_x[dim-1] < -m_overlap-m_epsilon || local_y[dim-1] < -m_overlap-m_epsilon)
412 if (nx*ny > m_max_normal_product + m_epsilon)
416 auto& intersection = m_edge_intersections[m_number_of_edge_intersections++];
417 intersection = { {{edgex, edgey}}, {{local_x, local_y}} };
419 catch(
const Dune::FMatrixError&) {
426 template<
typename Coordinate>
427 template<
typename Corners,
typename Normals>
429 ::projectionFeasible(
const Coordinate& x,
const Coordinate& nx,
const Coordinate& px,
const Corners& corners,
const Normals& normals)
const 431 using namespace ProjectionImplementation;
434 if (!
inside(px, m_epsilon))
438 if (px[dim-1] < -m_overlap-m_epsilon)
445 const auto d = xmy * n;
446 if (d < -m_overlap-m_epsilon)
450 if (nx * n > m_max_normal_product + m_epsilon)
457 template<
typename Coordinate>
458 template<
typename Corners,
typename Normals>
460 ::project(
const std::tuple<Corners&, Corners&>& corners,
const std::tuple<Normals&, Normals&>& normals)
462 doProjection(corners, normals);
463 doInverseProjection(corners, normals);
464 doEdgeIntersection(corners, normals);
Normals::value_type interpolate_unit_normals(const Coordinate &x, const Normals &normals)
Definition: projection_impl.hh:89
bool inside(const Coordinate &x, const Field &epsilon)
Definition: projection_impl.hh:109
Projection of a line (triangle) on another line (triangle).
Definition: projection.hh:18
Coordinate::field_type Field
Scalar type.
Definition: projection.hh:54
Corners::value_type interpolate(const Coordinate &x, const Corners &corners)
Definition: projection_impl.hh:68
std::pair< unsigned, unsigned > edgeToCorners(unsigned edge)
Definition: projection_impl.hh:41
Projection(const Field overlap=Field(0), const Field max_normal_product=Field(-0.1))
Definition: projection_impl.hh:128
void project(const std::tuple< Corners &, Corners &> &corners, const std::tuple< Normals &, Normals &> &normals)
Do the actual projection.
Definition: projection_impl.hh:460
void epsilon(const Field epsilon)
Set epsilon used for floating-point comparisons.
Definition: projection_impl.hh:138
Coordinate corner(unsigned c)
Definition: projection_impl.hh:22
Definition: gridglue.hh:30