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

Remove lex_faces from matrix-free kernels #13909

Merged
merged 1 commit into from
Jun 10, 2022
Merged

Remove lex_faces from matrix-free kernels #13909

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
20 changes: 7 additions & 13 deletions include/deal.II/matrix_free/evaluation_kernels.h
View file
Open in desktop
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 @@ -3841,17 +3838,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 @@ -3863,8 +3858,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
Open in desktop
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 @@ -447,8 +441,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 @@ -459,21 +452,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 stride = Utilities::pow(n_rows, face_direction);
Expand Down Expand Up @@ -534,55 +521,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 @@ -734,8 +707,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 @@ -899,15 +871,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 @@ -2505,8 +2474,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 @@ -2618,8 +2586,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 @@ -2631,12 +2598,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