Skip to content

Commit

Permalink
Use ndarray instead of c-array
Browse files Browse the repository at this point in the history 
Update comments and documentation
  • Loading branch information
@NiklasWik
NiklasWik committed May 2, 2022
1 parent 0661dbb commit a3b2bce
Show file tree
Hide file tree
Showing 4 changed files with 48 additions and 36 deletions.
11 changes: 5 additions & 6 deletions include/deal.II/matrix_free/evaluation_kernels.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@

#include <deal.II/base/config.h>

#include <deal.II/base/ndarray.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/vectorization.h>

Expand Down Expand Up @@ -3158,9 +3159,8 @@ namespace internal
(std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
n_dofs_tangent;

const unsigned int component_table[3][3] = {{1, 2, 0},
{2, 0, 1},
{0, 1, 2}};
static constexpr dealii::ndarray<unsigned int, 3, 3> component_table = {
{{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
const unsigned int component =
(dim == 2 && normal_dir == 0 && face_direction == 1) ?
0 :
Expand Down Expand Up @@ -3335,9 +3335,8 @@ namespace internal
(std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
n_dofs_tangent;

const unsigned int component_table[3][3] = {{1, 2, 0},
{2, 0, 1},
{0, 1, 2}};
static constexpr dealii::ndarray<unsigned int, 3, 3> component_table = {
{{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
const unsigned int component =
(dim == 2 && normal_dir == 0 && face_direction == 1) ?
0 :
Expand Down
63 changes: 33 additions & 30 deletions include/deal.II/matrix_free/fe_evaluation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5725,10 +5725,10 @@ inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_value(
const unsigned int q_point) const
{
// Check if Piola transform is required
if (this->data->element_type ==
internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
{
// Piola transform is required
# ifdef DEBUG
Assert(this->values_quad_initialized == true,
internal::ExcAccessToUninitializedField());
Expand All @@ -5741,23 +5741,22 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_value(
const std::size_t nqp = this->n_quadrature_points;
Tensor<1, dim, VectorizedArrayType> value_out;

// Cartesian cell
if (!is_face &&
this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
// Cartesian cell
const Tensor<2, dim, dealii::VectorizedArray<Number>> jac =
this->jacobian[1];
const VectorizedArrayType inv_det = determinant(this->jacobian[0]);

// J * u * det(J^-1)
for (unsigned int comp = 0; comp < n_components; ++comp)
value_out[comp] = this->values_quad[comp * nqp + q_point] *
jac[comp][comp] *
inv_det; // / this->jacobian[0][comp][comp];
jac[comp][comp] * inv_det;
}

// Affine or general cell
else
{
// Affine or general cell
const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
(this->cell_type > internal::MatrixFreeFunctions::affine) ?
this->jacobian[q_point] :
Expand Down Expand Up @@ -5786,6 +5785,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_value(
}
else
{
// No Piola needed
return BaseClass::get_value(q_point);
}
}
Expand All @@ -5795,10 +5795,10 @@ inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
get_gradient(const unsigned int q_point) const
{
// Check if Piola transform is required
if (this->data->element_type ==
internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
{
// Piola transform is required
# ifdef DEBUG
Assert(this->gradients_quad_initialized == true,
internal::ExcAccessToUninitializedField());
Expand All @@ -5811,24 +5811,25 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
const std::size_t nqp = this->n_quadrature_points;
Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_out;

// Cartesian cell
if (!is_face &&
this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
// Cartesian cell
const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
const VectorizedArrayType inv_det = determinant(inv_t_jac);

// J * grad_quad * J^-1 * det(J^-1)
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int comp = 0; comp < n_components; ++comp)
grad_out[comp][d] =
this->gradients_quad[(comp * dim + d) * nqp + q_point] *
inv_t_jac[d][d] * jac[comp][comp] * inv_det;
}
// Affine cell
else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
{
// Affine cell
const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
Expand All @@ -5853,9 +5854,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
grad_out[comp][d] = tmp;
}
}
// General cell TODO
else
{
// General cell
// Here we need the jacobian gradient and not the inverse which is
// stored in this->jacobian_gradients
AssertThrow(false, ExcNotImplemented());
Expand Down Expand Up @@ -5890,40 +5891,41 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
if (this->data->element_type ==
internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
{
// Affine cell
if (this->cell_type <= internal::MatrixFreeFunctions::affine)
{
// Affine cell
// Derivatives are reordered for faces. Need to take this into account
const VectorizedArrayType inv_det =
(is_face && dim == 2 && this->get_face_no() < 2) ?
-determinant(this->jacobian[0]) :
determinant(this->jacobian[0]);

// div * det(J^-1)
divergence = this->gradients_quad[q_point] * inv_det;
for (unsigned int d = 1; d < dim; ++d)
divergence +=
this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
}
// General cell TODO
else
{
// General cell
Assert(false, ExcNotImplemented());
}
}
else
{
// Cartesian cell
if (!is_face &&
this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
// Cartesian cell
divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
for (unsigned int d = 1; d < dim; ++d)
divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
this->jacobian[0][d][d];
}
// cell with general/constant Jacobian
else
{
// cell with general/constant Jacobian
const Tensor<2, dim, VectorizedArrayType> &jac =
this->cell_type == internal::MatrixFreeFunctions::general ?
this->jacobian[q_point] :
Expand Down Expand Up @@ -6039,13 +6041,11 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
submit_value(const Tensor<1, dim, VectorizedArrayType> val_in,
const unsigned int q_point)
{
// Check if Piola transform is required
if (this->data->element_type ==
internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
{
// Piola transform is required
AssertIndexRange(q_point, this->n_quadrature_points);

// This is not needed, but might be good to check anyway?
Assert(this->J_value != nullptr,
internal::ExcMatrixFreeAccessToUninitializedMappingField(
"update_value"));
Expand All @@ -6066,9 +6066,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
this->values_quad[comp * nqp + q_point] =
val_in[comp] * weight * jac[comp][comp];
}
// Affine or general cell
else
{
// Affine or general cell
const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
(this->cell_type > internal::MatrixFreeFunctions::affine) ?
this->jacobian[q_point] :
Expand All @@ -6090,7 +6090,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
-determinant(inv_t_jac) :
determinant(inv_t_jac)));

// J^T * u * w
// J^T * u * factor
for (unsigned int comp = 0; comp < n_components; ++comp)
{
this->values_quad[comp * nqp + q_point] =
Expand All @@ -6103,6 +6103,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
}
else
{
// No Piola transform
BaseClass::submit_value(val_in, q_point);
}
}
Expand All @@ -6113,10 +6114,11 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
submit_gradient(const Tensor<2, dim, VectorizedArrayType> grad_in,
const unsigned int q_point)
{
// Check if Piola transform is required
if (this->data->element_type ==
internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
{
// Piola transform is required

# ifdef DEBUG
Assert(this->is_reinitialized, ExcNotInitialized());
# endif
Expand All @@ -6132,10 +6134,10 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
# endif

const std::size_t nqp = this->n_quadrature_points;
// Cartesian cell
if (!is_face &&
this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
// Cartesian cell
const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
Expand All @@ -6145,9 +6147,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
this->gradients_quad[(comp * dim + d) * nqp + q_point] =
grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight;
}
// Affine cell
else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
{
// Affine cell
const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
Expand All @@ -6173,9 +6175,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp;
}
}
// General cell TODO
else
{
// General cell
AssertThrow(false, ExcNotImplemented());
}
}
Expand All @@ -6194,10 +6196,11 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
const unsigned int q_point)
{
// Check if Piola transform is required
if (this->data->element_type ==
internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
{
// Piola transform is required

# ifdef DEBUG
Assert(this->is_reinitialized, ExcNotInitialized());
# endif
Expand All @@ -6213,10 +6216,10 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
# endif

const std::size_t nqp = this->n_quadrature_points;
// Cartesian cell
if (!is_face &&
this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
// Cartesian cell
const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
Expand All @@ -6226,9 +6229,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
this->gradients_quad[(comp * dim + d) * nqp + q_point] =
grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight;
}
// Affine cell
else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
{
// Affine cell
const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
this->jacobian[0];
const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
Expand All @@ -6254,9 +6257,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp;
}
}
// General cell TODO
else
{
// General cell
AssertThrow(false, ExcNotImplemented());
}
}
Expand Down Expand Up @@ -6292,9 +6295,10 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
if (this->data->element_type ==
internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
{
// Affine cell
if (this->cell_type <= internal::MatrixFreeFunctions::affine)
{
// Affine cell

// Derivatives are reordered for faces. Need to take this into account
// and 1/inv_det != J_value for faces
const VectorizedArrayType fac =
Expand All @@ -6317,9 +6321,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
}
}
}
// General cell TODO
else
{
// General cell
AssertThrow(false, ExcNotImplemented());
}
}
Expand Down Expand Up @@ -6373,7 +6377,6 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
const SymmetricTensor<2, dim, VectorizedArrayType> sym_grad,
const unsigned int q_point)
{
// TODO
AssertThrow(
this->data->element_type !=
internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas,
Expand Down
2 changes: 2 additions & 0 deletions include/deal.II/matrix_free/fe_evaluation_data.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -743,6 +743,8 @@ class FEEvaluationData
* 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]. If the Jacobian also is stored,
* the components are instead reordered in the same way.
* Filled from MappingInfoStorage.jacobians in
* include/deal.II/matrix_free/mapping_info.templates.h
*/
const Tensor<2, dim, Number> *jacobian;

Expand Down
8 changes: 8 additions & 0 deletions include/deal.II/matrix_free/mapping_info_storage.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -227,6 +227,14 @@ namespace internal
* Contains two fields for access from both sides for interior faces,
* but the default case (cell integrals or boundary integrals) only
* fills the zeroth component and ignores the first one.
*
* If the cell is Cartesian/affine then the Jacobian is stored at index 1
* of the AlignedVector. For faces on hypercube elements, the derivatives
* are reorder s.t the derivative orthogonal to the face is stored last,
* i.e for dim = 3 and 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]. If the Jacobian also is stored, the components are
* instead reordered in the same way.
*/
std::array<AlignedVector<Tensor<2, spacedim, Number>>, 2> jacobians;

Expand Down

0 comments on commit a3b2bce

Please sign in to comment.