Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Avoid using Tensor::begin_raw() et al. #14647

Merged
merged 1 commit into from
Jan 9, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
11 changes: 5 additions & 6 deletions include/deal.II/base/bounding_box.h
Original file line number Diff line number Diff line change
Expand Up @@ -432,12 +432,11 @@ inline BoundingBox<spacedim, Number>::BoundingBox(const Container &points)
{
auto &min = boundary_points.first;
auto &max = boundary_points.second;
std::fill(min.begin_raw(),
min.end_raw(),
std::numeric_limits<Number>::infinity());
std::fill(max.begin_raw(),
max.end_raw(),
-std::numeric_limits<Number>::infinity());
for (unsigned int d = 0; d < spacedim; ++d)
{
min[d] = std::numeric_limits<Number>::infinity();
max[d] = -std::numeric_limits<Number>::infinity();
}

for (const Point<spacedim, Number> &point : points)
for (unsigned int d = 0; d < spacedim; ++d)
Expand Down
17 changes: 7 additions & 10 deletions source/base/quadrature_lib.cc
Original file line number Diff line number Diff line change
Expand Up @@ -615,13 +615,11 @@ template <>
unsigned int
QGaussOneOverR<2>::quad_size(const Point<2> &singularity, const unsigned int n)
{
const double eps = 1e-8;
const bool on_edge =
std::any_of(singularity.begin_raw(),
singularity.end_raw(),
[eps](double coord) {
return std::abs(coord) < eps || std::abs(coord - 1.) < eps;
});
const double eps = 1e-8;
bool on_edge = false;
for (unsigned int d = 0; d < 2; ++d)
on_edge = on_edge || (std::abs(singularity[d]) < eps ||
std::abs(singularity[d] - 1.0) < eps);
const bool on_vertex =
on_edge &&
std::abs((singularity - Point<2>(.5, .5)).norm_square() - .5) < eps;
Expand Down Expand Up @@ -2022,9 +2020,8 @@ QWitherdenVincentSimplex<dim>::QWitherdenVincentSimplex(
const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
this->weights.emplace_back(volume * b_weights[permutation_n]);
Point<dim> c_point;
std::copy(b_point.begin(),
b_point.begin() + dim,
c_point.begin_raw());
for (int d = 0; d < dim; ++d)
c_point[d] = b_point[d];
this->quadrature_points.emplace_back(c_point);
}
}
Expand Down
17 changes: 14 additions & 3 deletions source/base/tensor.cc
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
//
// ---------------------------------------------------------------------

#include <deal.II/base/array_view.h>
#include <deal.II/base/tensor.h>

#include <deal.II/lac/exceptions.h>
Expand Down Expand Up @@ -45,22 +46,32 @@ namespace
const types::blas_int lwork = 5 * dim;
std::array<Number, lwork> work;
types::blas_int info;
constexpr std::size_t size =
Tensor<2, dim, Number>::n_independent_components;
std::array<Number, size> A_array;
A_in_VT_out.unroll(A_array.begin(), A_array.end());
std::array<Number, size> U_array;
U.unroll(U_array.begin(), U_array.end());
gesvd(&LAPACKSupport::O, // replace VT in place
&LAPACKSupport::A,
&N,
&N,
A_in_VT_out.begin_raw(),
A_array.data(),
&N,
S.data(),
A_in_VT_out.begin_raw(),
A_array.data(),
&N,
U.begin_raw(),
U_array.data(),
&N,
work.data(),
&lwork,
&info);
Assert(info == 0, LAPACKSupport::ExcErrorCode("gesvd", info));
Assert(S.back() / S.front() > 1.e-10, LACExceptions::ExcSingular());

A_in_VT_out =
Tensor<2, dim, Number>(make_array_view(A_array.begin(), A_array.end()));
U = Tensor<2, dim, Number>(make_array_view(U_array.begin(), U_array.end()));
}
} // namespace

Expand Down
7 changes: 2 additions & 5 deletions source/fe/fe_simplex_p.cc
Original file line number Diff line number Diff line change
Expand Up @@ -95,11 +95,8 @@ namespace
// centroid and only the centroid
if (degree == 0)
{
Point<dim> centroid;
std::fill(centroid.begin_raw(),
centroid.end_raw(),
1.0 / double(dim + 1));
unit_points.emplace_back(centroid);
unit_points.emplace_back(
ReferenceCells::get_simplex<dim>().template barycenter<dim>());
return unit_points;
}

Expand Down
8 changes: 4 additions & 4 deletions source/fe/fe_simplex_p_bubbles.cc
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,8 @@ namespace FE_P_BubblesImplementation
FE_SimplexP<dim> fe_p(degree);
std::vector<Point<dim>> points = fe_p.get_unit_support_points();

Point<dim> centroid;
std::fill(centroid.begin_raw(), centroid.end_raw(), 1.0 / double(dim + 1));
const Point<dim> centroid =
fe_p.reference_cell().template barycenter<dim>();

switch (dim)
{
Expand Down Expand Up @@ -117,8 +117,8 @@ namespace FE_P_BubblesImplementation
BarycentricPolynomials<dim>
get_basis(const unsigned int degree)
{
Point<dim> centroid;
std::fill(centroid.begin_raw(), centroid.end_raw(), 1.0 / double(dim + 1));
const Point<dim> centroid =
ReferenceCells::get_simplex<dim>().template barycenter<dim>();

auto M = [](const unsigned int d) {
return BarycentricPolynomial<dim, double>::monomial(d);
Expand Down
6 changes: 4 additions & 2 deletions tests/base/quadrature_point_data_04.cc
Original file line number Diff line number Diff line change
Expand Up @@ -66,13 +66,15 @@ struct Mat1 : MaterialBase
pack_values(std::vector<double> &values) const final
{
AssertDimension(values.size(), number_of_values());
std::copy(pt.begin_raw(), pt.end_raw(), values.begin());
for (unsigned int d = 0; d < 2; ++d)
values[d] = pt[d];
}
virtual void
unpack_values(const std::vector<double> &values) final
{
AssertDimension(values.size(), number_of_values());
std::copy(values.cbegin(), values.cend(), pt.begin_raw());
for (unsigned int d = 0; d < 2; ++d)
pt[d] = values[d];
}
};

Expand Down
21 changes: 13 additions & 8 deletions tests/symengine/symengine_tensor_operations_04.cc
Original file line number Diff line number Diff line change
Expand Up @@ -104,10 +104,12 @@ test_tensor()
using Tensor_t = Tensor<rank, dim, double>;

Tensor_t t_a, t_b;
for (auto it = t_a.begin_raw(); it != t_a.end_raw(); ++it)
*it = 1.0;
for (auto it = t_b.begin_raw(); it != t_b.end_raw(); ++it)
*it = 2.0;
for (unsigned int i = 0; i < Tensor_t::n_independent_components; ++i)
{
const auto index = t_a.unrolled_to_component_index(i);
t_a[index] = 1.0;
t_b[index] = 2.0;
}

const Tensor_SD_number_t symb_t_a =
SD::make_tensor_of_symbols<rank, dim>("a");
Expand All @@ -128,10 +130,13 @@ test_symmetric_tensor()
using Tensor_t = SymmetricTensor<rank, dim, double>;

Tensor_t t_a, t_b;
for (auto it = t_a.begin_raw(); it != t_a.end_raw(); ++it)
*it = 1.0;
for (auto it = t_b.begin_raw(); it != t_b.end_raw(); ++it)
*it = 2.0;
for (unsigned int i = 0; i < Tensor_t::n_independent_components; ++i)
{
const auto index = t_a.unrolled_to_component_index(i);
t_a[index] = 1.0;
t_b[index] = 2.0;
}


const Tensor_SD_number_t symb_t_a =
SD::make_symmetric_tensor_of_symbols<rank, dim>("a");
Expand Down