From b584c58d9558ec8fba128a098404499825388d14 Mon Sep 17 00:00:00 2001 From: David Wells Date: Mon, 12 Jun 2017 11:13:04 -0400 Subject: [PATCH] Improve vertex caching in a few places. Calling cell->vertex() requires several index lookups before indexing into the global vertex array so this is relatively expensive. This patch lowers, in a TransfiniteInterpolationManifold benchmark, the number of vertex lookups from about 361 million to 229 million and lowers the total wall time by about 5%. --- source/grid/manifold_lib.cc | 20 ++++++++++++++------ source/grid/tria_accessor.cc | 2 +- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 491052f53d..5566da7ea8 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -694,12 +694,14 @@ namespace // https://en.wikipedia.org/wiki/Transfinite_interpolation // S(u,v) = (1-v)c_1(u)+v c_3(u) + (1-u)c_2(v) + u c_4(v) - // [(1-u)(1-v)P_0 + u(1-v) P_1 + (1-u)v P_2 + uv P_3] + const std::array, GeometryInfo<2>::vertices_per_cell> vertices + {{cell.vertex(0), cell.vertex(1), cell.vertex(2), cell.vertex(3)}}; 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) * - cell.vertex(v); + vertices[v]; else { // We subtract the contribution of the vertices (second line in formula). @@ -733,8 +735,8 @@ namespace } else { - points[0] = cell.vertex(GeometryInfo<2>::line_to_cell_vertices(line,0)); - points[1] = cell.vertex(GeometryInfo<2>::line_to_cell_vertices(line,1)); + points[0] = vertices[GeometryInfo<2>::line_to_cell_vertices(line,0)]; + points[1] = vertices[GeometryInfo<2>::line_to_cell_vertices(line,1)]; weights[0] = 1. - line_point; weights[1] = line_point; new_point += my_weight * @@ -744,7 +746,7 @@ namespace // subtract contribution from the vertices (second line in formula) for (unsigned int v=0; v::vertices_per_cell; ++v) - new_point += weights_vertices[v] * cell.vertex(v); + new_point += weights_vertices[v] * vertices[v]; } return new_point; @@ -1009,16 +1011,22 @@ TransfiniteInterpolationManifold if (&cell->get_manifold() != this) continue; + std::array, GeometryInfo::vertices_per_cell> vertices; + for (unsigned int vertex_n = 0; vertex_n < GeometryInfo::vertices_per_cell; ++vertex_n) + { + vertices[vertex_n] = cell->vertex(vertex_n); + } + // cheap check: if any of the points is not inside a circle around the // center of the loop, we can skip the expensive part below (this assumes // that the manifold does not deform the grid too much) Point center; for (unsigned int v=0; v::vertices_per_cell; ++v) - center += cell->vertex(v); + center += vertices[v]; center *= 1./GeometryInfo::vertices_per_cell; double radius_square = 0.; for (unsigned int v=0; v::vertices_per_cell; ++v) - radius_square = std::max(radius_square, (center-cell->vertex(v)).norm_square()); + radius_square = std::max(radius_square, (center-vertices[v]).norm_square()); bool inside_circle = true; for (unsigned int i=0; i radius_square * 1.5) diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 207fc909e1..e23ed92602 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1349,7 +1349,7 @@ TriaAccessor // b = vertex * Kb Tensor<1,spacedim> b = point; for (unsigned int v=0; v::vertices_per_cell; ++v) - b -= this->vertex(v) * TransformR2UAffine::Kb[v]; + b -= vertices[v] * TransformR2UAffine::Kb[v]; DerivativeForm<1,spacedim,structdim> A_inv = A.covariant_form().transpose(); return Point(apply_transformation(A_inv, b)); -- 2.39.5