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.

Sign up for GitHub

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

Jump to bottom

Make all template specialization of FPE compile for float #14155

Closed
wants to merge 1 commit into from
Closed

Make all template specialization of FPE compile for float #14155

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
14 changes: 9 additions & 5 deletions include/deal.II/base/derivative_form.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,12 +108,14 @@ class DerivativeForm
* Number>. In particular, if order == 1 and the derivative is the Jacobian of
* $\mathbf F(\mathbf x)$, then Tensor[i] = $\nabla F_i(\mathbf x)$.
*/
operator Tensor<order + 1, dim, Number>() const;
template <typename Number2>
operator Tensor<order + 1, dim, Number2>() const;

/**
* Converts a DerivativeForm<1, dim, 1, Number> to Tensor<1, dim, Number>.
*/
operator Tensor<1, dim, Number>() const;
template <typename Number2>
operator Tensor<1, dim, Number2>() const;

/**
* Return the transpose of a rectangular DerivativeForm,
Expand Down Expand Up @@ -273,8 +275,9 @@ DerivativeForm<order, dim, spacedim, Number>::operator[](


template <int order, int dim, int spacedim, typename Number>
template <typename Number2>
inline DerivativeForm<order, dim, spacedim, Number>::
operator Tensor<1, dim, Number>() const
operator Tensor<1, dim, Number2>() const
{
Assert((1 == spacedim) && (order == 1),
ExcMessage("Only allowed for spacedim==1."));
Expand All @@ -285,12 +288,13 @@ operator Tensor<1, dim, Number>() const


template <int order, int dim, int spacedim, typename Number>
template <typename Number2>
inline DerivativeForm<order, dim, spacedim, Number>::
operator Tensor<order + 1, dim, Number>() const
operator Tensor<order + 1, dim, Number2>() const
Comment on lines +291 to +293
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This doesn't sound quite right to me. Shouldn't the operator just return a Tensor representation of what it stores internally, and then if someone wants this in another type they can later cast from that tensor to one with a different scalar?

Note that that's basically what you do in the implementation, which now looks as follows:

template <int order, int dim, int spacedim, typename Number>
template <typename Number2>
inline DerivativeForm<order, dim, spacedim, Number>::
operator Tensor<order + 1, dim, Number>() const
{
  Assert((dim == spacedim), ExcMessage("Only allowed when dim==spacedim."));

  Tensor<order + 1, dim, Number> t;       <<<<<<<<<<<<<<<<<< old tensor type
  if (dim == spacedim)
    for (unsigned int j = 0; j < dim; ++j)
      t[j] = (*this)[j];
  return t;                                        <<<<<<<<<<<< cast to Tensor<....,Number2> happens here
}

My preference would be to make the place that calls this function do the cast at the calling site, not internally here.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tend to agree with @bangerth here, since this only converts the end result, I would prefer doing that explicitly in FEPointEvaluation (I was confused by the abbreviation FPE, which I read as floating point exception).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That the tensor used the old number- type was not indended:

Tensor<order + 1, dim, Number> t;       <<<<<<<<<<<<<<<<<< old tensor type

I have adjusted that.

Since the operator() is for casting, the casting is already done explicitly done in FEPointEvaluation:

dealii/include/deal.II/matrix_free/fe_point_evaluation.h

Lines 950 to 954 in 1195dbb

gradients[i + j] = static_cast<gradient_type>(
apply_transformation(mapping_info->get_mapping_data()
.inverse_jacobians[i + j]
.transpose(),
unit_gradients[i + j]));

The alternative would to cast from DerivativeForm to typename internal::FEPointEvaluation::EvaluatorTypeTraits<dim, n_components, double>::gradient_type (please note that one has to use double instead of Number in this long expression) and after that cast it to a float version (or do an assignment, which seems to broke in Tensor).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think what I meant to say is this: I consider the cast operator as saying "give me a tensor representation of this object", the emphasis being on "representation of this object" rather than "give me". The emphasis here is on the object. As a consequence, a DerivativeForm<...,double> should return a Tensor<...,double>.

In contrast, your implementation is for the question "I want to use this object as an X", in which the emphasis is on want to use -- a caller side consideration. I would prefer an explicit cast to Tensor<...,Number2> in such situations.

{
Assert((dim == spacedim), ExcMessage("Only allowed when dim==spacedim."));

Tensor<order + 1, dim, Number> t;
Tensor<order + 1, dim, Number2> t;

if (dim == spacedim)
for (unsigned int j = 0; j < dim; ++j)
Expand Down
18 changes: 6 additions & 12 deletions include/deal.II/matrix_free/fe_point_evaluation.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -947,15 +947,11 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
Number>::set_gradient(val_and_grad.second,
j,
unit_gradients[i + j]);
gradients[i + j] =
static_cast<typename internal::FEPointEvaluation::
EvaluatorTypeTraits<dim,
n_components,
Number>::gradient_type>(
apply_transformation(mapping_info->get_mapping_data()
.inverse_jacobians[i + j]
.transpose(),
unit_gradients[i + j]));
gradients[i + j] = static_cast<gradient_type>(
apply_transformation(mapping_info->get_mapping_data()
.inverse_jacobians[i + j]
.transpose(),
unit_gradients[i + j]));
}
}
}
Expand Down Expand Up @@ -1092,9 +1088,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
{
gradients[i + j] =
static_cast<typename internal::FEPointEvaluation::
EvaluatorTypeTraits<dim, n_components, Number>::
gradient_type>(apply_transformation(
static_cast<gradient_type>(apply_transformation(
mapping_info->get_mapping_data().inverse_jacobians[i + j],
gradients[i + j]));
internal::FEPointEvaluation::
Expand Down
76 changes: 76 additions & 0 deletions tests/matrix_free/point_evaluation_18.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
// ---------------------------------------------------------------------
//
// Copyright (C) 2022 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


// Check that all template specializations of FEPointEvaluation are
// compiling.

#include <deal.II/matrix_free/fe_point_evaluation.h>

#include "../tests.h"

template <int n_components, int dim, int spacedim, typename Number>
void
test()
{
return; // nothing to do, since we are only interested if the code
// compiles

std::unique_ptr<Mapping<dim, spacedim>> mapping;
std::unique_ptr<FiniteElement<dim, spacedim>> fe;

FEPointEvaluation<n_components, dim, spacedim, Number> fpe(
*mapping, *fe, UpdateFlags::update_default);

Triangulation<dim, spacedim> tria;

fpe.reinit(tria.begin(), ArrayView<const Point<dim>>());

fpe.evaluate(ArrayView<const Number>(), EvaluationFlags::values);

fpe.integrate(ArrayView<Number>(), EvaluationFlags::values);
}

int
main()
{
initlog();

test<1, 1, 1, double>();
test<2, 1, 1, double>();

test<1, 2, 2, double>();
test<2, 2, 2, double>();
test<3, 2, 2, double>();

test<1, 3, 3, double>();
test<2, 3, 3, double>();
test<3, 3, 3, double>();
test<4, 3, 3, double>();

test<1, 1, 1, float>();
test<2, 1, 1, float>();

test<1, 2, 2, float>();
test<2, 2, 2, float>();
test<3, 2, 2, float>();

test<1, 3, 3, float>();
test<2, 3, 3, float>();
test<3, 3, 3, float>();
test<4, 3, 3, float>();

deallog << "OK!" << std::endl;
}
2 changes: 2 additions & 0 deletions tests/matrix_free/point_evaluation_18.output
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@

DEAL::OK!