Skip to content

Commit

Permalink
Merge pull request #13909 from peterrum/lex_faces_remove
Browse files Browse the repository at this point in the history 
Remove lex_faces from matrix-free kernels
  • Loading branch information
@kronbichler
kronbichler committed Jun 10, 2022
2 parents 2d52c32 + 5a99157 commit 3abca26
Show file tree
Hide file tree
Showing 2 changed files with 48 additions and 90 deletions.
20 changes: 7 additions & 13 deletions include/deal.II/matrix_free/evaluation_kernels.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3536,7 +3536,7 @@ namespace internal
};


template <int dim, int fe_degree, typename Number, bool lex_faces = false>
template <int dim, int fe_degree, typename Number>
struct FEFaceNormalEvaluationImpl
{
template <bool do_evaluate, bool add_into_output>
Expand Down Expand Up @@ -3637,20 +3637,17 @@ namespace internal
evalf.template apply_face<face_direction,
do_evaluate,
add_into_output,
2,
lex_faces>(input, output);
2>(input, output);
else if (flag & EvaluationFlags::gradients)
evalf.template apply_face<face_direction,
do_evaluate,
add_into_output,
1,
lex_faces>(input, output);
1>(input, output);
else
evalf.template apply_face<face_direction,
do_evaluate,
add_into_output,
0,
lex_faces>(input, output);
0>(input, output);
input += in_stride;
output += out_stride;
}
Expand Down Expand Up @@ -3849,17 +3846,15 @@ namespace internal
evalf0.template apply_face<face_direction,
do_evaluate,
add_into_output,
max_derivative,
lex_faces>(input, output);
max_derivative>(input, output);
// stride to next component
input += (face_direction == 0) ? in_stride_after_normal : in_stride;
output += (face_direction == 0) ? out_stride_after_normal : out_stride;

evalf1.template apply_face<face_direction,
do_evaluate,
add_into_output,
max_derivative,
lex_faces>(input, output);
max_derivative>(input, output);

if (dim == 3)
{
Expand All @@ -3871,8 +3866,7 @@ namespace internal
evalf2.template apply_face<face_direction,
do_evaluate,
add_into_output,
max_derivative,
lex_faces>(input, output);
max_derivative>(input, output);
}
}
};
Expand Down
118 changes: 41 additions & 77 deletions include/deal.II/matrix_free/tensor_product_kernels.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -334,20 +334,14 @@ namespace internal
* derivatives, 2 second derivates. Note that all the
* derivatives access the data in @p shape_values passed to
* the constructor of the class
* @tparam lex_faces Sets how the evaluation points on the faces should be
* sorted: lexicographically or right-hand-system number
* (special treatment of orientation 1 in 3D). Per default
* right-hand-system number is enabled, which is only
* working for dimensions up to 3.
*
* @param in address of the input data vector
* @param out address of the output data vector
*/
template <int face_direction,
bool contract_onto_face,
bool add,
int max_derivative,
bool lex_faces = false>
int max_derivative>
void
apply_face(const Number *DEAL_II_RESTRICT in,
Number *DEAL_II_RESTRICT out) const;
Expand Down Expand Up @@ -448,8 +442,7 @@ namespace internal
template <int face_direction,
bool contract_onto_face,
bool add,
int max_derivative,
bool lex_faces>
int max_derivative>
inline void
EvaluatorTensorProduct<evaluate_general,
dim,
Expand All @@ -460,21 +453,15 @@ namespace internal
Number *DEAL_II_RESTRICT
out) const
{
Assert(dim > 0 && (lex_faces || dim < 4),
ExcMessage("Only dim=1,2,3 supported"));
Assert(dim > 0, ExcMessage("Only dim=1,2,3 supported"));
static_assert(max_derivative >= 0 && max_derivative < 3,
"Only derivative orders 0-2 implemented");
Assert(shape_values != nullptr,
ExcMessage(
"The given array shape_values must not be the null pointer."));

constexpr int n_blocks1 =
lex_faces ? dealii::Utilities::pow<unsigned int>(n_rows, face_direction) :
(dim > 1 ? n_rows : 1);
constexpr int n_blocks2 =
lex_faces ? dealii::Utilities::pow<unsigned int>(
n_rows, std::max(dim - face_direction - 1, 0)) :
(dim > 2 ? n_rows : 1);
constexpr int n_blocks1 = (dim > 1 ? n_rows : 1);
constexpr int n_blocks2 = (dim > 2 ? n_rows : 1);

AssertIndexRange(face_direction, dim);
constexpr int in_stride = Utilities::pow(n_rows, face_direction);
Expand Down Expand Up @@ -536,55 +523,41 @@ namespace internal
}
}

if (lex_faces)
// increment: in regular case, just go to the next point in
// x-direction. If we are at the end of one chunk in x-dir, need
// to jump over to the next layer in z-direction
switch (face_direction)
{
++out;
++in;
case 0:
in += contract_onto_face ? n_rows : 1;
out += contract_onto_face ? 1 : n_rows;
break;
case 1:
++in;
++out;
// faces 2 and 3 in 3D use local coordinate system zx, which
// is the other way around compared to the tensor
// product. Need to take that into account.
if (dim == 3)
{
if (contract_onto_face)
out += n_rows - 1;
else
in += n_rows - 1;
}
break;
case 2:
++in;
++out;
break;
default:
Assert(false, ExcNotImplemented());
}
else
// increment: in regular case, just go to the next point in
// x-direction. If we are at the end of one chunk in x-dir, need
// to jump over to the next layer in z-direction
switch (face_direction)
{
case 0:
in += contract_onto_face ? n_rows : 1;
out += contract_onto_face ? 1 : n_rows;
break;
case 1:
++in;
++out;
// faces 2 and 3 in 3D use local coordinate system zx, which
// is the other way around compared to the tensor
// product. Need to take that into account.
if (dim == 3)
{
if (contract_onto_face)
out += n_rows - 1;
else
in += n_rows - 1;
}
break;
case 2:
++in;
++out;
break;
default:
Assert(false, ExcNotImplemented());
}
}
if (lex_faces)
{
if (contract_onto_face)
in += (dealii::Utilities::pow(n_rows, face_direction + 1) -
n_blocks1);
else
out += (dealii::Utilities::pow(n_rows, face_direction + 1) -
n_blocks1);
}
else if (face_direction == 1 && dim == 3)

// adjust for local coordinate system zx
if (face_direction == 1 && dim == 3)
{
// adjust for local coordinate system zx
if (contract_onto_face)
{
in += n_rows * (n_rows - 1);
Expand Down Expand Up @@ -736,8 +709,7 @@ namespace internal
template <int face_direction,
bool contract_onto_face,
bool add,
int max_derivative,
bool lex_faces = false>
int max_derivative>
void
apply_face(const Number *DEAL_II_RESTRICT in,
Number *DEAL_II_RESTRICT out) const;
Expand Down Expand Up @@ -901,15 +873,12 @@ namespace internal
template <int face_direction,
bool contract_onto_face,
bool add,
int max_derivative,
bool lex_faces>
int max_derivative>
inline void
EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::
apply_face(const Number *DEAL_II_RESTRICT in,
Number *DEAL_II_RESTRICT out) const
{
static_assert(lex_faces == false, "Not implemented yet.");

Assert(shape_values != nullptr,
ExcMessage(
"The given array shape_data must not be the null pointer!"));
Expand Down Expand Up @@ -2511,8 +2480,7 @@ namespace internal
template <int face_direction,
bool contract_onto_face,
bool add,
int max_derivative,
bool lex_faces = false>
int max_derivative>
void
apply_face(const Number *DEAL_II_RESTRICT in,
Number *DEAL_II_RESTRICT out) const;
Expand Down Expand Up @@ -2625,8 +2593,7 @@ namespace internal
template <int face_direction,
bool contract_onto_face,
bool add,
int max_derivative,
bool lex_faces>
int max_derivative>
inline void
EvaluatorTensorProductAnisotropic<
evaluate_raviart_thomas,
Expand All @@ -2638,12 +2605,9 @@ namespace internal
Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
Number *DEAL_II_RESTRICT out) const
{
Assert(dim > 1 && (lex_faces || dim < 4),
ExcMessage("Only dim=2,3 supported"));
Assert(dim > 1 && dim < 4, ExcMessage("Only dim=2,3 supported"));
static_assert(max_derivative >= 0 && max_derivative < 3,
"Only derivative orders 0-2 implemented");
static_assert(!lex_faces,
"lex_faces = True is not implemented for Raviart-Thomas");
Assert(shape_values != nullptr,
ExcMessage(
"The given array shape_values must not be the null pointer."));
Expand Down

0 comments on commit 3abca26

Please sign in to comment.