Skip to content

Commit

Permalink
Avoid using Tensor::begin_raw() et al.
Browse files Browse the repository at this point in the history
These are deprecated and in several cases we have better choices anyway.
  • Loading branch information
drwells committed Jan 8, 2023
1 parent 4ca9bfe commit 7ad0689
Show file tree
Hide file tree
Showing 7 changed files with 49 additions and 38 deletions.
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

0 comments on commit 7ad0689

Please sign in to comment.