From 147f527d2d2999245fce70e14e2cce40c07968bb Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 10 Nov 2017 09:35:50 +0100 Subject: [PATCH] Further speedup of transfinite interpolation. Try to use an initial guess for pull_back that is assuming cube-like geometries. Skip some indirections inside compute_transfinite_interpolation. --- include/deal.II/grid/manifold_lib.h | 12 ++- source/grid/manifold_lib.cc | 155 +++++++++++++++++++--------- 2 files changed, 119 insertions(+), 48 deletions(-) diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index 0cf05a34b5..fa58b747e4 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -685,6 +685,10 @@ public: * this class is that the input triangulation is uniformly refined and the * manifold is later attached to the same triangulation. * + * Whenever the assignment of manifold ids changes on the level of the + * triangulation which this class was initialized with, initialize() must be + * called again to update the manifold ids connected to the coarse cells. + * * @note The triangulation used to construct the manifold must not be * destroyed during the usage of this object. */ @@ -764,6 +768,11 @@ private: /** * Pull back operation into the unit coordinates on the given coarse cell. * + * This method is currently based on a Newton-like iteration to find the + * point in the origin. One may speed up the iteration by providing a good + * initial guess as the third argument. If no better point is known, use + * cell->real_to_unit_cell_affine_approximation(p) + * * @note This internal function is currently not compatible with the * ChartManifold::pull_back() function because the given class represents an * atlas of charts, not a single chart. Thus, the pull_back() operation is @@ -774,7 +783,8 @@ private: */ Point pull_back(const typename Triangulation::cell_iterator &cell, - const Point &p) const; + const Point &p, + const Point &initial_guess) const; /** * Push forward operation. diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index ebac5f2764..f7b6eea6cc 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -929,8 +929,10 @@ namespace const Point<2> &chart_point, const bool cell_is_flat) { + const unsigned int dim = AccessorType::dimension; const unsigned int spacedim = AccessorType::space_dimension; const types::manifold_id my_manifold_id = cell.manifold_id(); + const Triangulation &tria = cell.get_triangulation(); // formula see wikipedia // https://en.wikipedia.org/wiki/Transfinite_interpolation @@ -939,11 +941,20 @@ namespace const std::array, 4> vertices {{cell.vertex(0), cell.vertex(1), cell.vertex(2), cell.vertex(3)}}; + // store the components of the linear shape functions because we need them + // repeatedly + double linear_shapes[4]; + for (unsigned int d=0; d<2; ++d) + { + linear_shapes[2*d] = 1.-chart_point[d]; + linear_shapes[2*d+1] = chart_point[d]; + } + Point new_point; if (cell_is_flat) - for (unsigned int v=0; v::vertices_per_cell; ++v) - new_point += GeometryInfo<2>::d_linear_shape_function(chart_point, v) * - vertices[v]; + for (unsigned int i1=0, v=0; i1<2; ++i1) + for (unsigned int i0=0; i0<2; ++i0, ++v) + new_point += linear_shapes[2+i1] * linear_shapes[i0] * vertices[v]; else { // We subtract the contribution of the vertices (second line in formula). @@ -953,8 +964,9 @@ namespace // in 2D but it becomes clear in 3D where we avoid looking at the faces' // orientation and other complications). double weights_vertices[GeometryInfo<2>::vertices_per_cell]; - for (unsigned int v=0; v::vertices_per_cell; ++v) - weights_vertices[v] = -GeometryInfo<2>::d_linear_shape_function(chart_point, v); + for (unsigned int i1=0, v=0; i1<2; ++i1) + for (unsigned int i0=0; i0<2; ++i0, ++v) + weights_vertices[v] = -linear_shapes[2+i1] * linear_shapes[i0]; // add the contribution from the lines around the cell (first line in // formula) @@ -971,8 +983,9 @@ namespace // Same manifold or invalid id which will go back to the same class -> // adds to the vertices - if (cell.line(line)->manifold_id() == my_manifold_id || - cell.line(line)->manifold_id() == numbers::invalid_manifold_id) + const types::manifold_id line_manifold_id = cell.line(line)->manifold_id(); + if (line_manifold_id == my_manifold_id || + line_manifold_id == numbers::invalid_manifold_id) { weights_vertices[GeometryInfo<2>::line_to_cell_vertices(line,0)] += my_weight * (1.-line_point); @@ -986,8 +999,8 @@ namespace weights[0] = 1. - line_point; weights[1] = line_point; new_point += my_weight * - cell.line(line)->get_manifold().get_new_point(points_view, - weights_view); + tria.get_manifold(line_manifold_id).get_new_point(points_view, + weights_view); } } @@ -1037,6 +1050,7 @@ namespace const unsigned int dim = AccessorType::dimension; const unsigned int spacedim = AccessorType::space_dimension; const types::manifold_id my_manifold_id = cell.manifold_id(); + const Triangulation &tria = cell.get_triangulation(); // Same approach as in 2D, but adding the faces, subtracting the edges, and // adding the vertices @@ -1047,18 +1061,34 @@ namespace cell.vertex(4), cell.vertex(5), cell.vertex(6), cell.vertex(7) } }; - Point new_point; + // store the components of the linear shape functions because we need them + // repeatedly. we allow for 10 such shape functions to wrap around the + // first four once again for easier face access. + double linear_shapes[10]; + for (unsigned int d=0; d<3; ++d) + { + linear_shapes[2*d] = 1.-chart_point[d]; + linear_shapes[2*d+1] = chart_point[d]; + } + + Point new_point; if (cell_is_flat) - for (unsigned int v=0; v::vertices_per_cell; ++v) - new_point += GeometryInfo<3>::d_linear_shape_function(chart_point, v) * - vertices[v]; + for (unsigned int i2=0, v=0; i2<2; ++i2) + for (unsigned int i1=0; i1<2; ++i1) + for (unsigned int i0=0; i0<2; ++i0, ++v) + new_point += (linear_shapes[4+i2] * linear_shapes[2+i1]) * + linear_shapes[i0] * vertices[v]; else { // identify the weights for the vertices and lines to be accumulated double weights_vertices[GeometryInfo<3>::vertices_per_cell]; - for (unsigned int v=0; v::vertices_per_cell; ++v) - weights_vertices[v] = GeometryInfo<3>::d_linear_shape_function(chart_point, v); + for (unsigned int i2=0, v=0; i2<2; ++i2) + for (unsigned int i1=0; i1<2; ++i1) + for (unsigned int i0=0; i0<2; ++i0, ++v) + weights_vertices[v] = (linear_shapes[4+i2] * linear_shapes[2+i1]) * + linear_shapes[i0]; + double weights_lines[GeometryInfo<3>::lines_per_cell]; for (unsigned int line=0; line::lines_per_cell; ++line) weights_lines[line] = 0; @@ -1070,39 +1100,50 @@ namespace const auto weights_view = make_array_view(weights.begin(), weights.end()); const auto points_view = make_array_view(points.begin(), points.end()); + // wrap linear shape functions around for access in face loop + for (unsigned int d=6; d<10; ++d) + linear_shapes[d] = linear_shapes[d-6]; + for (unsigned int face=0; face::faces_per_cell; ++face) { - Point<2> quad_point(chart_point[(face/2+1)%3], chart_point[(face/2+2)%3]); - const double my_weight = face%2 ? chart_point[face/2] : 1-chart_point[face/2]; + const double my_weight = linear_shapes[face]; + const unsigned int face_even = face - face%2; if (std::abs(my_weight) < 1e-13) continue; // same manifold or invalid id which will go back to the same class // -> face will interpolate from the surrounding lines and vertices - if (cell.face(face)->manifold_id() == my_manifold_id || - cell.face(face)->manifold_id() == numbers::invalid_manifold_id) + const types::manifold_id face_manifold_id = cell.face(face)->manifold_id(); + if (face_manifold_id == my_manifold_id || + face_manifold_id == numbers::invalid_manifold_id) { for (unsigned int line=0; line::lines_per_cell; ++line) { - const double line_weight = line%2 ? quad_point[line/2] : 1-quad_point[line/2]; + const double line_weight = linear_shapes[face_even+2+line]; weights_lines[face_to_cell_lines_3d[face][line]] += my_weight * line_weight; } - for (unsigned int v=0; v::vertices_per_cell; ++v) - weights_vertices[face_to_cell_vertices_3d[face][v]] - -= GeometryInfo<2>::d_linear_shape_function(quad_point, v) * my_weight; + weights_vertices[face_to_cell_vertices_3d[face][0]] -= + linear_shapes[face_even+2]*(linear_shapes[face_even+4]*my_weight); + weights_vertices[face_to_cell_vertices_3d[face][1]] -= + linear_shapes[face_even+3]*(linear_shapes[face_even+4]*my_weight); + weights_vertices[face_to_cell_vertices_3d[face][2]] -= + linear_shapes[face_even+2]*(linear_shapes[face_even+5]*my_weight); + weights_vertices[face_to_cell_vertices_3d[face][3]] -= + linear_shapes[face_even+3]*(linear_shapes[face_even+5]*my_weight); } else { for (unsigned int v=0; v::vertices_per_cell; ++v) - { - points[v] = vertices[face_to_cell_vertices_3d[face][v]]; - weights[v] = GeometryInfo<2>::d_linear_shape_function(quad_point, v); - } + points[v] = vertices[face_to_cell_vertices_3d[face][v]]; + weights[0] = linear_shapes[face_even+2]*linear_shapes[face_even+4]; + weights[1] = linear_shapes[face_even+3]*linear_shapes[face_even+4]; + weights[2] = linear_shapes[face_even+2]*linear_shapes[face_even+5]; + weights[3] = linear_shapes[face_even+3]*linear_shapes[face_even+5]; new_point += my_weight * - cell.face(face)->get_manifold().get_new_point(points_view, - weights_view); + tria.get_manifold(face_manifold_id).get_new_point(points_view, + weights_view); } } @@ -1114,23 +1155,20 @@ namespace const double line_point = (line < 8 ? chart_point[1-(line%4)/2] : chart_point[2]); double my_weight = 0.; if (line < 8) - { - const unsigned int subline = line%4; - my_weight = subline % 2 ? chart_point[subline/2] : 1-chart_point[subline/2]; - my_weight *= line/4 ? chart_point[2] : (1-chart_point[2]); - } + my_weight = linear_shapes[line%4] * linear_shapes[4+line/4]; else { - Point<2> xy(chart_point[0], chart_point[1]); - my_weight = GeometryInfo<2>::d_linear_shape_function(xy, line-8); + const unsigned int subline = line-8; + my_weight = linear_shapes[subline%2] * linear_shapes[2+subline/2]; } my_weight -= weights_lines[line]; if (std::abs(my_weight) < 1e-13) continue; - if (cell.line(line)->manifold_id() == my_manifold_id || - cell.line(line)->manifold_id() == numbers::invalid_manifold_id) + const types::manifold_id line_manifold_id = cell.line(line)->manifold_id(); + if (line_manifold_id == my_manifold_id || + line_manifold_id == numbers::invalid_manifold_id) { weights_vertices[GeometryInfo<3>::line_to_cell_vertices(line,0)] -= my_weight * (1.-line_point); @@ -1144,8 +1182,8 @@ namespace weights[0] = 1. - line_point; weights[1] = line_point; new_point -= my_weight * - cell.line(line)->get_manifold().get_new_point(points_view_line, - weights_view_line); + tria.get_manifold(line_manifold_id).get_new_point(points_view_line, + weights_view_line); } } @@ -1211,15 +1249,15 @@ template Point TransfiniteInterpolationManifold ::pull_back(const typename Triangulation::cell_iterator &cell, - const Point &point) const + const Point &point, + const Point &initial_guess) const { Point outside; for (unsigned int d=0; d chart_point = - GeometryInfo::project_to_unit_cell(cell->real_to_unit_cell_affine_approximation(point)); + // project the user-given input to unit cell + Point chart_point = GeometryInfo::project_to_unit_cell(initial_guess); // run quasi-Newton iteration with a combination of finite differences for // the exact Jacobian and "Broyden's good method". As opposed to the various @@ -1228,7 +1266,7 @@ TransfiniteInterpolationManifold // the implementation of the compute_chart_points method. Tensor<1,spacedim> residual = point - compute_transfinite_interpolation(*cell, chart_point, coarse_cell_is_flat[cell->index()]); - const double tolerance = 1e-21 * cell->diameter() * cell->diameter(); + const double tolerance = 1e-21 * Utilities::fixed_power<2>(cell->diameter()); double residual_norm_square = residual.norm_square(); DerivativeForm<1,dim,spacedim> inv_grad; for (unsigned int i=0; i<100; ++i) @@ -1426,7 +1464,29 @@ TransfiniteInterpolationManifold bool inside_unit_cell = true; for (unsigned int i=0; i p3 = chart_points[1] + + Point(chart_points[2]-chart_points[0]); + chart_points[i] = pull_back(cell, surrounding_points[i], p3); + } + // 8 points usually form either a cube or a rectangle with vertices + // and line mid points. assume a cube here which gives us some new + // initial guesses + else if (surrounding_points.size() == 8 && i > 4) + { + const Point guess = chart_points[i-4] + + Point(chart_points[4] - chart_points[0]); + chart_points[i] = pull_back(cell, surrounding_points[i], guess); + } + else + chart_points[i] = pull_back(cell, surrounding_points[i], + cell->real_to_unit_cell_affine_approximation(surrounding_points[i])); // Tolerance 1e-6 chosen that the method also works with // SphericalManifold @@ -1464,7 +1524,8 @@ TransfiniteInterpolationManifold for (unsigned int i=0; i " - << pull_back(cell, surrounding_points[i]) + << pull_back(cell, surrounding_points[i], + cell->real_to_unit_cell_affine_approximation(surrounding_points[i])) << std::endl; } } -- 2.39.5