Skip to content

Commit

Permalink
Avoid computing inverse jac in piola
Browse files Browse the repository at this point in the history 
Updated documentation and also storing jacobian for Cartesian/affine faces
  • Loading branch information
@NiklasWik
NiklasWik committed Apr 28, 2022
1 parent da0dd32 commit 2dacda6
Show file tree
Hide file tree
Showing 3 changed files with 61 additions and 21 deletions.
28 changes: 16 additions & 12 deletions include/deal.II/matrix_free/fe_evaluation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5762,7 +5762,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_value(
this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
const Tensor<2, dim, dealii::VectorizedArray<Number>> jac =
invert(this->jacobian[0]);
this->jacobian[1];
const VectorizedArrayType inv_det = determinant(this->jacobian[0]);

for (unsigned int comp = 0; comp < n_components; ++comp)
Expand All @@ -5778,7 +5778,10 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_value(
(this->cell_type > internal::MatrixFreeFunctions::affine) ?
this->jacobian[q_point] :
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
const Tensor<2, dim, VectorizedArrayType> &t_jac =
(this->cell_type > internal::MatrixFreeFunctions::affine) ?
invert(inv_t_jac) :
this->jacobian[1];

// Derivatives are reordered for faces. Need to take this into account
const std::uint8_t fn = (is_face) ? this->get_face_no() : 0;
Expand Down Expand Up @@ -5831,7 +5834,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
{
const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
const Tensor<2, dim, VectorizedArrayType> &t_jac = this->jacobian[1];
const VectorizedArrayType inv_det = determinant(inv_t_jac);

for (unsigned int d = 0; d < dim; ++d)
Expand All @@ -5845,7 +5848,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
{
const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
const Tensor<2, dim, VectorizedArrayType> &t_jac = this->jacobian[1];

// Derivatives are reordered for faces. Need to take this into account
const std::uint8_t fn = (is_face) ? this->get_face_no() : 0;
Expand Down Expand Up @@ -6076,7 +6079,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
const Tensor<2, dim, dealii::VectorizedArray<Number>> t_jac =
invert(this->jacobian[0]);
this->jacobian[1];
const VectorizedArrayType weight = this->quadrature_weights[q_point];

for (unsigned int comp = 0; comp < n_components; ++comp)
Expand All @@ -6090,7 +6093,10 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
(this->cell_type > internal::MatrixFreeFunctions::affine) ?
this->jacobian[q_point] :
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
const Tensor<2, dim, VectorizedArrayType> &t_jac =
(this->cell_type > internal::MatrixFreeFunctions::affine) ?
invert(inv_t_jac) :
this->jacobian[1];

// Derivatives are reordered for faces. Need to take this into account
const std::uint8_t fn = (is_face) ? this->get_face_no() : 0;
Expand Down Expand Up @@ -6154,8 +6160,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
{
const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &t_jac =
invert(this->jacobian[0]);
const Tensor<2, dim, VectorizedArrayType> &t_jac = this->jacobian[1];
const VectorizedArrayType weight = this->quadrature_weights[q_point];
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int comp = 0; comp < n_components; ++comp)
Expand All @@ -6167,7 +6172,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
{
const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
const Tensor<2, dim, VectorizedArrayType> &t_jac = this->jacobian[1];

// Derivatives are reordered for faces. Need to take this into account
const std::uint8_t fn = (is_face) ? this->get_face_no() : 0;
Expand Down Expand Up @@ -6239,8 +6244,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
{
const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &t_jac =
invert(this->jacobian[0]);
const Tensor<2, dim, VectorizedArrayType> &t_jac = this->jacobian[1];
const VectorizedArrayType weight = this->quadrature_weights[q_point];
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int comp = 0; comp < n_components; ++comp)
Expand All @@ -6252,7 +6256,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
{
const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
const Tensor<2, dim, VectorizedArrayType> &t_jac = this->jacobian[1];

// Derivatives are reordered for faces. Need to take this into account
const std::uint8_t fn = (is_face) ? this->get_face_no() : 0;
Expand Down
20 changes: 13 additions & 7 deletions include/deal.II/matrix_free/fe_evaluation_data.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -735,8 +735,13 @@ class FEEvaluationData
const Point<dim, Number> *quadrature_points;

/**
* A pointer to the Jacobian information of the present cell. Only set to a
* useful value if on a non-Cartesian cell.
* A pointer to the inverse transpose Jacobian information of the present
* cell. Only set to a useful value if on a non-Cartesian cell. If the cell is
* Cartesian/affine then the transpose Jacobian is stored at index 1. For
* faces, the derivatives are reorder s.t the normal derivative is stored
* last, i.e for face_no = 0 or 1, the derivatives are ordered as
* [dy, dz, dx], face_no = 2 or 3: [dz, dx, dy], and face_no = 5 or 6:
* [dx, dy, dz].
*/
const Tensor<2, dim, Number> *jacobian;

Expand Down Expand Up @@ -1279,16 +1284,17 @@ FEEvaluationData<dim, Number, is_face>::quadrature_point(
Assert(this->jacobian != nullptr, ExcNotInitialized());
Point<dim, Number> point = this->quadrature_points[0];

const Tensor<2, dim, Number> &jac = this->jacobian[1];
const Tensor<2, dim, Number> &t_jac = this->jacobian[1];
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
for (unsigned int d = 0; d < dim; ++d)
point[d] += jac[d][d] * static_cast<typename Number::value_type>(
this->descriptor->quadrature.point(q)[d]);
point[d] += t_jac[d][d] * static_cast<typename Number::value_type>(
this->descriptor->quadrature.point(q)[d]);
else
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int e = 0; e < dim; ++e)
point[d] += jac[d][e] * static_cast<typename Number::value_type>(
this->descriptor->quadrature.point(q)[e]);
point[d] +=
t_jac[e][d] * static_cast<typename Number::value_type>(
this->descriptor->quadrature.point(q)[e]);
return point;
}
else
Expand Down
34 changes: 32 additions & 2 deletions include/deal.II/matrix_free/mapping_info.templates.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1214,7 +1214,7 @@ namespace internal
store_vectorized_array(
jac[d][e],
vv,
my_data.jacobians[0][idx + 1][d][e]);
my_data.jacobians[0][idx + 1][e][d]);
}
else
{
Expand Down Expand Up @@ -2177,6 +2177,21 @@ namespace internal
vv,
my_data.jacobians[0][offset + q][d][e]);
}
// Also store the jacobian
// Flip columns instead so that it matches if you call
// invert(jacobians[q]) for general faces.
if (face_type[face] <= affine)
for (unsigned int e = 0; e < dim; ++e)
{
const unsigned int ee =
ExtractFaceHelper::reorder_face_derivative_indices<
dim>(interior_face_no, e);
for (unsigned int d = 0; d < dim; ++d)
store_vectorized_array(
jac[d][ee],
vv,
my_data.jacobians[0][offset + q + 1][e][d]);
}

if (update_flags_faces & update_jacobian_grads)
{
Expand Down Expand Up @@ -2281,6 +2296,21 @@ namespace internal
vv,
my_data.jacobians[1][offset + q][d][e]);
}
// Also store the jacobian
// Flip columns instead so that it matches if you call
// invert(jacobians[q]) for general faces.
if (face_type[face] <= affine)
for (unsigned int e = 0; e < dim; ++e)
{
const unsigned int ee = ExtractFaceHelper::
reorder_face_derivative_indices<dim>(
exterior_face_no, e);
for (unsigned int d = 0; d < dim; ++d)
store_vectorized_array(
jac[d][ee],
vv,
my_data.jacobians[1][offset + q + 1][e][d]);
}

if (update_flags_faces & update_jacobian_grads)
{
Expand Down Expand Up @@ -2762,7 +2792,7 @@ namespace internal
max_size =
std::max(max_size,
my_data.data_index_offsets[face] +
(face_type[face] <= affine ? 1 : n_q_points));
(face_type[face] <= affine ? 2 : n_q_points));
}

const UpdateFlags update_flags_common =
Expand Down

0 comments on commit 2dacda6

Please sign in to comment.