From c4d593de59c5f364f6d52d6b9a3ad2b14e95fce7 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 1 Apr 2024 21:20:23 +0200 Subject: [PATCH] MappingCartesian: Some performance improvements --- include/deal.II/fe/mapping_cartesian.h | 16 ++- source/fe/mapping_cartesian.cc | 185 ++++++++++++------------- 2 files changed, 107 insertions(+), 94 deletions(-) diff --git a/include/deal.II/fe/mapping_cartesian.h b/include/deal.II/fe/mapping_cartesian.h index 390091b5a3..73f1f2006f 100644 --- a/include/deal.II/fe/mapping_cartesian.h +++ b/include/deal.II/fe/mapping_cartesian.h @@ -109,6 +109,13 @@ public: const typename Triangulation::cell_iterator &cell, const Point &p) const override; + // for documentation, see the Mapping base class + virtual void + transform_points_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const ArrayView> &real_points, + const ArrayView> &unit_points) const override; + /** * @} */ @@ -223,6 +230,13 @@ public: */ mutable Tensor<1, dim> cell_extents; + /** + * Reciprocal of the extents of the last cell we have seen in the + * coordinate directions, i.e., hx, + * hy, hz. + */ + mutable Tensor<1, dim> inverse_cell_extents; + /** * The volume element */ @@ -352,7 +366,7 @@ private: /** * Transform quadrature points in InternalData to real space by scaling unit - * coordinates with cell_extends in each direction. + * coordinates with cell_extents in each direction. * * Called from the various maybe_update_*_quadrature_points functions. */ diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc index 85be4007ee..947e473df6 100644 --- a/source/fe/mapping_cartesian.cc +++ b/source/fe/mapping_cartesian.cc @@ -91,6 +91,7 @@ template MappingCartesian::InternalData::InternalData( const Quadrature &q) : cell_extents(numbers::signaling_nan>()) + , inverse_cell_extents(numbers::signaling_nan>()) , volume_element(numbers::signaling_nan()) , quadrature_points(q.get_points()) {} @@ -103,6 +104,7 @@ MappingCartesian::InternalData::memory_consumption() const { return (Mapping::InternalDataBase::memory_consumption() + MemoryConsumption::memory_consumption(cell_extents) + + MemoryConsumption::memory_consumption(inverse_cell_extents) + MemoryConsumption::memory_consumption(volume_element) + MemoryConsumption::memory_consumption(quadrature_points)); } @@ -232,29 +234,18 @@ MappingCartesian::update_cell_extents( const CellSimilarity::Similarity cell_similarity, const InternalData &data) const { - // Compute start point and sizes - // along axes. Strange vertex - // numbering makes this complicated - // again. + // Compute start point and sizes along axes. The vertices to be looked at + // are 1, 2, 4 compared to the base vertex 0. if (cell_similarity != CellSimilarity::translation) { - const Point &start = cell->vertex(0); - switch (dim) + const Point start = cell->vertex(0); + for (unsigned int d = 0; d < dim; ++d) { - case 1: - data.cell_extents[0] = cell->vertex(1)[0] - start[0]; - break; - case 2: - data.cell_extents[0] = cell->vertex(1)[0] - start[0]; - data.cell_extents[1] = cell->vertex(2)[1] - start[1]; - break; - case 3: - data.cell_extents[0] = cell->vertex(1)[0] - start[0]; - data.cell_extents[1] = cell->vertex(2)[1] - start[1]; - data.cell_extents[2] = cell->vertex(4)[2] - start[2]; - break; - default: - DEAL_II_NOT_IMPLEMENTED(); + const double cell_extent_d = cell->vertex(1 << d)[d] - start[d]; + data.cell_extents[d] = cell_extent_d; + Assert(cell_extent_d != 0., + ExcMessage("Cell does not appear to be Cartesian!")); + data.inverse_cell_extents[d] = 1. / cell_extent_d; } } } @@ -345,7 +336,7 @@ MappingCartesian::transform_quadrature_points( { // let @p{start} be the origin of a local coordinate system. it is chosen // as the (lower) left vertex - const Point &start = cell->vertex(0); + const Point start = cell->vertex(0); for (unsigned int i = 0; i < quadrature_points.size(); ++i) { @@ -483,7 +474,8 @@ MappingCartesian::maybe_update_inverse_jacobians( { output_data.inverse_jacobians[i] = Tensor<2, dim>(); for (unsigned int j = 0; j < dim; ++j) - output_data.inverse_jacobians[i][j][j] = 1. / data.cell_extents[j]; + output_data.inverse_jacobians[i][j][j] = + data.inverse_cell_extents[j]; } } @@ -721,7 +713,7 @@ MappingCartesian::fill_fe_immersed_surface_values( const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i); for (unsigned int d = 0; d < dim; ++d) { - normal[d] = ref_space_normal[d] / data.cell_extents[d]; + normal[d] = ref_space_normal[d] * data.inverse_cell_extents[d]; } normal /= normal.norm(); output_data.normal_vectors[i] = normal; @@ -738,7 +730,8 @@ MappingCartesian::fill_fe_immersed_surface_values( for (unsigned int d = 0; d < dim; ++d) { det_jacobian *= data.cell_extents[d]; - invJTxNormal[d] = ref_space_normal[d] / data.cell_extents[d]; + invJTxNormal[d] = + ref_space_normal[d] * data.inverse_cell_extents[d]; } output_data.JxW_values[i] = det_jacobian * invJTxNormal.norm() * quadrature.weight(i); @@ -775,7 +768,7 @@ MappingCartesian::transform( for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d = 0; d < dim; ++d) - output[i][d] = input[i][d] / data.cell_extents[d]; + output[i][d] = input[i][d] * data.inverse_cell_extents[d]; return; } @@ -836,7 +829,8 @@ MappingCartesian::transform( for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) - output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2]; + output[i][d1][d2] = + input[i][d1][d2] * data.inverse_cell_extents[d2]; return; } @@ -862,8 +856,9 @@ MappingCartesian::transform( for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) - output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2] / - data.cell_extents[d1]; + output[i][d1][d2] = input[i][d1][d2] * + data.inverse_cell_extents[d2] * + data.inverse_cell_extents[d1]; return; } @@ -876,8 +871,8 @@ MappingCartesian::transform( for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) - output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] / - data.cell_extents[d1]; + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * + data.inverse_cell_extents[d1]; return; } @@ -910,8 +905,9 @@ MappingCartesian::transform( for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) - output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] / - data.cell_extents[d1] / data.volume_element; + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * + data.inverse_cell_extents[d1] / + data.volume_element; return; } @@ -946,7 +942,8 @@ MappingCartesian::transform( for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) - output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2]; + output[i][d1][d2] = + input[i][d1][d2] * data.inverse_cell_extents[d2]; return; } @@ -972,8 +969,9 @@ MappingCartesian::transform( for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) - output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2] / - data.cell_extents[d1]; + output[i][d1][d2] = input[i][d1][d2] * + data.inverse_cell_extents[d2] * + data.inverse_cell_extents[d1]; return; } @@ -986,8 +984,8 @@ MappingCartesian::transform( for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) - output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] / - data.cell_extents[d1]; + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * + data.inverse_cell_extents[d1]; return; } @@ -1020,8 +1018,9 @@ MappingCartesian::transform( for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) - output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] / - data.cell_extents[d1] / data.volume_element; + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * + data.inverse_cell_extents[d1] / + data.volume_element; return; } @@ -1058,9 +1057,9 @@ MappingCartesian::transform( for (unsigned int j = 0; j < spacedim; ++j) for (unsigned int k = 0; k < spacedim; ++k) { - output[q][i][j][k] = input[q][i][j][k] / - data.cell_extents[j] / - data.cell_extents[k]; + output[q][i][j][k] = input[q][i][j][k] * + data.inverse_cell_extents[j] * + data.inverse_cell_extents[k]; } return; } @@ -1100,9 +1099,10 @@ MappingCartesian::transform( for (unsigned int j = 0; j < spacedim; ++j) for (unsigned int k = 0; k < spacedim; ++k) { - output[q][i][j][k] = - input[q][i][j][k] * data.cell_extents[i] / - data.cell_extents[j] / data.cell_extents[k]; + output[q][i][j][k] = input[q][i][j][k] * + data.cell_extents[i] * + data.inverse_cell_extents[j] * + data.inverse_cell_extents[k]; } return; } @@ -1118,9 +1118,10 @@ MappingCartesian::transform( for (unsigned int j = 0; j < spacedim; ++j) for (unsigned int k = 0; k < spacedim; ++k) { - output[q][i][j][k] = - input[q][i][j][k] / data.cell_extents[i] / - data.cell_extents[j] / data.cell_extents[k]; + output[q][i][j][k] = input[q][i][j][k] * + (data.inverse_cell_extents[i] * + data.inverse_cell_extents[j]) * + data.inverse_cell_extents[k]; } return; @@ -1144,9 +1145,10 @@ MappingCartesian::transform( for (unsigned int k = 0; k < spacedim; ++k) { output[q][i][j][k] = - input[q][i][j][k] * data.cell_extents[i] / - data.volume_element / data.cell_extents[j] / - data.cell_extents[k]; + input[q][i][j][k] * + (data.cell_extents[i] / data.volume_element * + data.inverse_cell_extents[j]) * + data.inverse_cell_extents[k]; } return; @@ -1166,32 +1168,16 @@ MappingCartesian::transform_unit_to_real_cell( const Point &p) const { Assert(is_cartesian(cell), ExcCellNotCartesian()); + if (dim != spacedim) + DEAL_II_NOT_IMPLEMENTED(); - Tensor<1, dim> length; - const Point start = cell->vertex(0); - switch (dim) - { - case 1: - length[0] = cell->vertex(1)[0] - start[0]; - break; - case 2: - length[0] = cell->vertex(1)[0] - start[0]; - length[1] = cell->vertex(2)[1] - start[1]; - break; - case 3: - length[0] = cell->vertex(1)[0] - start[0]; - length[1] = cell->vertex(2)[1] - start[1]; - length[2] = cell->vertex(4)[2] - start[2]; - break; - default: - DEAL_II_NOT_IMPLEMENTED(); - } + Point unit = cell->vertex(0); - Point p_real = cell->vertex(0); + // Go through vertices with numbers 1, 2, 4 for (unsigned int d = 0; d < dim; ++d) - p_real[d] += length[d] * p[d]; + unit[d] += (cell->vertex(1 << d)[d] - unit[d]) * p[d]; - return p_real; + return unit; } @@ -1206,32 +1192,45 @@ MappingCartesian::transform_real_to_unit_cell( if (dim != spacedim) DEAL_II_NOT_IMPLEMENTED(); - const Point &start = cell->vertex(0); - Point real = p; - real -= start; + const Point start = cell->vertex(0); + Point real = p; + + // Go through vertices with numbers 1, 2, 4 + for (unsigned int d = 0; d < dim; ++d) + real[d] = (real[d] - start[d]) / (cell->vertex(1 << d)[d] - start[d]); - switch (dim) - { - case 1: - real[0] /= cell->vertex(1)[0] - start[0]; - break; - case 2: - real[0] /= cell->vertex(1)[0] - start[0]; - real[1] /= cell->vertex(2)[1] - start[1]; - break; - case 3: - real[0] /= cell->vertex(1)[0] - start[0]; - real[1] /= cell->vertex(2)[1] - start[1]; - real[2] /= cell->vertex(4)[2] - start[2]; - break; - default: - DEAL_II_NOT_IMPLEMENTED(); - } return real; } +template +void +MappingCartesian::transform_points_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const ArrayView> &real_points, + const ArrayView> &unit_points) const +{ + Assert(is_cartesian(cell), ExcCellNotCartesian()); + AssertDimension(real_points.size(), unit_points.size()); + + if (dim != spacedim) + DEAL_II_NOT_IMPLEMENTED(); + + const Point start = cell->vertex(0); + + // Go through vertices with numbers 1, 2, 4 + std::array inverse_lengths; + for (unsigned int d = 0; d < dim; ++d) + inverse_lengths[d] = 1. / (cell->vertex(1 << d)[d] - start[d]); + + for (unsigned int i = 0; i < real_points.size(); ++i) + for (unsigned int d = 0; d < dim; ++d) + unit_points[i][d] = (real_points[i][d] - start[d]) * inverse_lengths[d]; +} + + + template std::unique_ptr> MappingCartesian::clone() const -- 2.39.5