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;
178 for (
unsigned i = 0; i < dim-1; ++i) {
179 directions[i] = target[i+1] - target[0];
180 scales[i] = directions[i].infinity_norm();
181 directions[i] /= scales[i];
184 for (
unsigned i = 0; i < dim-1; ++i) {
185 for (
unsigned j = 0; j < dim; ++j) {
186 m[j][i] = directions[i][j];
190 m_projection_valid =
true;
194 for (
unsigned i = 0; i < origin.size(); ++i) {
195 for (
unsigned j = 0; j < dim; ++j)
196 m[j][dim-1] = origin_normals[i][j];
198 const Coordinate rhs = origin[i] - target[0];
204 for (
unsigned j = 0; j < dim-1; ++j)
207 y[dim-1] *= Field(-1);
209 const bool feasible = projectionFeasible(origin[i], origin_normals[i], y, target, target_normals);
210 success.set(i, feasible);
212 catch (
const Dune::FMatrixError&) {
213 success.set(i,
false);
214 m_projection_valid =
false;
219 template<
typename Coordinate>
220 template<
typename Corners,
typename Normals>
222 Projection<Coordinate>
223 ::doInverseProjection(
const std::tuple<Corners&, Corners&>& corners,
const std::tuple<Normals&, Normals&>& normals)
235 using namespace ProjectionImplementation;
237 typedef Dune::FieldMatrix<Field, dim-1, dim-1> Matrix;
238 typedef Dune::FieldVector<Field, dim-1> Vector;
243 if (!m_projection_valid) {
244 get<1>(m_success).reset();
248 const auto& images = get<0>(m_images);
249 const auto& target_corners = get<1>(corners);
250 auto& preimages = get<1>(m_images);
251 auto& success = get<1>(m_success);
253 std::array<Coordinate, dim> v;
254 for (
unsigned i = 0; i < dim-1; ++i) {
260 for (
unsigned i = 0; i < dim-1; ++i) {
261 for (
unsigned j = 0; j < dim-1; ++j) {
266 for (
unsigned i = 0; i < dim; ++i) {
268 v[dim-1] = target_corners[i];
269 v[dim-1] -=
interpolate(images[0], target_corners);
272 for (
unsigned j = 0; j < dim-1; ++j)
273 rhs[j] = v[dim-1]*v[j];
276 for (
unsigned j = 0; j < dim-1; ++j)
277 preimages[i][j] = z[j];
281 preimages[i][dim-1] = (x - target_corners[i]) * get<1>(normals)[i];
284 const bool feasible = projectionFeasible(target_corners[i], get<1>(normals)[i], preimages[i], get<0>(corners), get<0>(normals));
285 success.set(i, feasible);
289 template<
typename Coordinate>
290 template<
typename Corners,
typename Normals>
292 Projection<Coordinate>
293 ::doEdgeIntersection(
const std::tuple<Corners&, Corners&>& corners,
const std::tuple<Normals&, Normals&>& normals)
295 using namespace ProjectionImplementation;
298 m_number_of_edge_intersections = 0;
309 if (!m_projection_valid || get<0>(m_success).all() || get<1>(m_success).all()) {
313 const auto& images = get<0>(m_images);
314 const auto& ys = get<1>(corners);
325 for (
unsigned edgex = 0; edgex < dim; ++edgex) {
330 if (get<0>(m_success)[i] && get<0>(m_success)[j])
335 const auto pxjpxi = pxj - pxi;
337 typedef Dune::FieldMatrix<Field, dim-1, dim-1> Matrix;
338 typedef Dune::FieldVector<Field, dim-1> Vector;
340 for (
unsigned edgey = 0; edgey < dim; ++edgey) {
345 if (get<1>(m_success)[k] && get<1>(m_success)[l])
348 const auto ykyl = ys[k] - ys[l];
349 const auto ykpxi = ys[k] - pxi;
354 for (
unsigned m = 0; m < dim-1; ++m) {
355 const auto ym1y0 = ys[m+1] - ys[0];
356 mat[m][0] = pxjpxi * ym1y0;
357 mat[m][1] = ykyl * ym1y0;
358 rhs[m] = ykpxi * ym1y0;
367 if (!isfinite(z[0]) || !isfinite(z[1]))
371 if (z[0] < m_epsilon || z[0] > Field(1) - m_epsilon
372 || z[1] < m_epsilon || z[1] > Field(1) - m_epsilon)
375 Coordinate local_x = corner<Coordinate, Field>(i);
376 local_x.axpy(z[0], corner<Coordinate, Field>(j) - corner<Coordinate, Field>(i));
377 Coordinate local_y = corner<Coordinate, Field>(k);
378 local_y.axpy(z[1], corner<Coordinate, Field>(l) - corner<Coordinate, Field>(k));
381 if (!
inside(local_x, m_epsilon) || !
inside(local_y, m_epsilon))
389 local_x[dim-1] = -(xy*nx);
390 local_y[dim-1] = xy*ny;
392 if (local_x[dim-1] < -m_overlap-m_epsilon || local_y[dim-1] < -m_overlap-m_epsilon)
396 if (nx*ny > m_max_normal_product + m_epsilon)
400 auto& intersection = m_edge_intersections[m_number_of_edge_intersections++];
401 intersection = { {{edgex, edgey}}, {{local_x, local_y}} };
403 catch(
const Dune::FMatrixError&) {
410 template<
typename Coordinate>
411 template<
typename Corners,
typename Normals>
412 bool Projection<Coordinate>
413 ::projectionFeasible(
const Coordinate& x,
const Coordinate& nx,
const Coordinate& px,
const Corners& corners,
const Normals& normals)
const
415 using namespace ProjectionImplementation;
418 if (!
inside(px, m_epsilon))
422 if (px[dim-1] < -m_overlap-m_epsilon)
429 const auto d = xmy * n;
430 if (d < -m_overlap-m_epsilon)
434 if (nx * n > m_max_normal_product + m_epsilon)
441 template<
typename Coordinate>
442 template<
typename Corners,
typename Normals>
444 ::project(
const std::tuple<Corners&, Corners&>& corners,
const std::tuple<Normals&, Normals&>& normals)
446 doProjection(corners, normals);
447 doInverseProjection(corners, normals);
448 doEdgeIntersection(corners, normals);
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
Projection(const Field overlap=Field(0), const Field max_normal_product=Field(-0.1))
Definition: projection_impl.hh:128
Definition: gridglue.hh:33
Corners::value_type interpolate(const Coordinate &x, const Corners &corners)
Definition: projection_impl.hh:68
Coordinate corner(unsigned c)
Definition: projection_impl.hh:22
void epsilon(const Field epsilon)
Set epsilon used for floating-point comparisons.
Definition: projection_impl.hh:138
Projection of a line (triangle) on another line (triangle).
Definition: projection.hh:18
void project(const std::tuple< Corners &, Corners & > &corners, const std::tuple< Normals &, Normals & > &normals)
Do the actual projection.
Definition: projection_impl.hh:444
Coordinate::field_type Field
Scalar type.
Definition: projection.hh:54
Normals::value_type interpolate_unit_normals(const Coordinate &x, const Normals &normals)
Definition: projection_impl.hh:89