Skip to content

Commit

Permalink
Fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@peterrum
peterrum committed Apr 22, 2022
1 parent 0258a8d commit be8f0d4
Show file tree
Hide file tree
Showing 7 changed files with 186 additions and 34 deletions.
105 changes: 85 additions & 20 deletions include/deal.II/matrix_free/evaluation_kernels.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2868,7 +2868,7 @@ namespace internal
{
Assert((fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.get_is_interior_face() == false) == false,
fe_eval.is_interior_face() == false) == false,
ExcNotImplemented());

const unsigned int face_no = fe_eval.get_face_no();
Expand Down Expand Up @@ -2944,15 +2944,20 @@ namespace internal

if (fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.get_is_interior_face() == false)
fe_eval.is_interior_face() == false)
{
for (unsigned int v = 0; v < Number::size(); ++v)
{
// the loop breaks once an invalid_unsigned_int is hit for
// all cases except the exterior faces in the ECL loop (where
// some faces might be at the boundaries but others not)
if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
continue;
{
for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
++i)
temp[i][v] = 0;
continue;
}

FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
template interpolate<true, false>(n_components,
Expand Down Expand Up @@ -3015,7 +3020,7 @@ namespace internal

if (fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.get_is_interior_face() == false)
fe_eval.is_interior_face() == false)
{
for (unsigned int v = 0; v < Number::size(); ++v)
{
Expand All @@ -3034,7 +3039,7 @@ namespace internal
&fe_eval.get_shape_info().face_orientations_quad(
fe_eval.get_face_orientation(v), 0),
false,
Utilities::pow(n_q_points_1d, dim - 1),
shape_info.n_q_points_face,
&temp[0][0],
fe_eval.begin_values(),
fe_eval.begin_gradients(),
Expand Down Expand Up @@ -3071,15 +3076,20 @@ namespace internal
Number * values_dofs,
FEEvaluationData<dim, Number, true> & fe_eval)
{
const auto & shape_info = fe_eval.get_shape_info();
const auto & shape_data = shape_info.data.front();
const unsigned int face_no = fe_eval.get_face_no();
const unsigned int face_orientation = fe_eval.get_face_orientation();
const auto &shape_info = fe_eval.get_shape_info();
const auto &shape_data = shape_info.data.front();

if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
{
const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
Assert((fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.is_interior_face() == false) == false,
ExcNotImplemented());

const unsigned int face_no = fe_eval.get_face_no();
const unsigned int face_orientation = fe_eval.get_face_orientation();
const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];

using Eval =
EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
Expand Down Expand Up @@ -3150,12 +3160,41 @@ namespace internal
Number *temp = fe_eval.get_scratch_data().begin();
Number *scratch_data = temp + 3 * n_components * dofs_per_face;

if (face_orientation)
if (fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.is_interior_face() == false)
{
for (unsigned int v = 0; v < Number::size(); ++v)
{
// the loop breaks once an invalid_unsigned_int is hit for
// all cases except the exterior faces in the ECL loop (where
// some faces might be at the boundaries but others not)
if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
continue;

if (fe_eval.get_face_orientation(v) != 0)
adjust_for_face_orientation_per_lane(
dim,
n_components,
v,
integration_flag,
&fe_eval.get_shape_info().face_orientations_quad(
fe_eval.get_face_orientation(v), 0),
true,
shape_info.n_q_points_face,
&temp[0][0],
fe_eval.begin_values(),
fe_eval.begin_gradients(),
fe_eval.begin_hessians());
}
}
else if (fe_eval.get_face_orientation() != 0)
adjust_for_face_orientation(
dim,
n_components,
integration_flag,
&fe_eval.get_shape_info().face_orientations_quad(face_orientation, 0),
&fe_eval.get_shape_info().face_orientations_quad(
fe_eval.get_face_orientation(), 0),
true,
shape_info.n_q_points_face,
temp,
Expand Down Expand Up @@ -3201,13 +3240,39 @@ namespace internal
scratch_data,
subface_index);

FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
template interpolate<false, false>(n_components,
integration_flag,
shape_info,
temp,
values_dofs,
face_no);
if (fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.is_interior_face() == false)
{
for (unsigned int v = 0; v < Number::size(); ++v)
{
// the loop breaks once an invalid_unsigned_int is hit for
// all cases except the exterior faces in the ECL loop (where
// some faces might be at the boundaries but others not)
if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
continue;

FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
template interpolate<false, false>(n_components,
integration_flag,
shape_info,
values_dofs,
scratch_data,
fe_eval.get_face_no(v));

for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
++i)
temp[i][v] = scratch_data[i][v];
}
}
else
FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
template interpolate<false, false>(n_components,
integration_flag,
shape_info,
temp,
values_dofs,
fe_eval.get_face_no());
return false;
}
};
Expand Down
21 changes: 19 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 @@ -3039,8 +3039,8 @@ namespace internal
for (unsigned int e = 0; e < dim; ++e)
{
const unsigned int ee = ExtractFaceHelper::
reorder_face_derivative_indices<dim>(face,
e);
reorder_face_derivative_indices<dim>(
fe_val_neigh.get_face_number(), e);
face_data_by_cells[my_q]
.jacobians[1][offset][d][e][v] =
inv_jac[d][ee];
Expand Down Expand Up @@ -3081,6 +3081,23 @@ namespace internal
inv_jac[d][ee];
}
}
if (is_local && (update_flags & update_jacobians))
for (unsigned int q = 0; q < fe_val.n_quadrature_points;
++q)
{
DerivativeForm<1, dim, dim> inv_jac =
fe_val_neigh.jacobian(q).covariant_form();
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int e = 0; e < dim; ++e)
{
const unsigned int ee = ExtractFaceHelper::
reorder_face_derivative_indices<dim>(
fe_val_neigh.get_face_number(), e);
face_data_by_cells[my_q]
.jacobians[1][offset + q][d][e][v] =
inv_jac[d][ee];
}
}
if (update_flags & update_jacobian_grads)
{
Assert(false, ExcNotImplemented());
Expand Down
15 changes: 13 additions & 2 deletions tests/matrix_free/ecl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,14 +51,25 @@ template <int dim,
typename Number = double,
typename VectorizedArrayType = VectorizedArray<Number>>
void
test(const unsigned int n_refinements = 1,
test(const unsigned int geometry = 0,
const unsigned int n_refinements = 1,
const bool print_vector = true,
const MPI_Comm comm = MPI_COMM_SELF)
{
using VectorType = LinearAlgebra::distributed::Vector<Number>;

parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
GridGenerator::hyper_cube(tria, 0.0, 1.0, true);

if (geometry == 0)
GridGenerator::hyper_cube(tria, 0.0, 1.0, true);
else if (geometry == 1)
GridGenerator::hyper_shell(tria, Point<dim>(), 0.5, 1.0);
else if (geometry == 2)
GridGenerator::hyper_ball(tria);
else
Assert(false, ExcNotImplemented());

tria.reset_all_manifolds();

if (false)
{
Expand Down
6 changes: 4 additions & 2 deletions tests/matrix_free/ecl_01.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,13 +34,15 @@


// Compare the evaluation of a SIP Laplace operator with face- and element-
// centric loops on a structured grid.
// centric loops on a hypercube, hypershell, and hyperball.

int
main(int argc, char **argv)
{
Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);

mpi_initlog();
test<2, 2, 3, double, VectorizedArray<double>>();
test<2, 2, 3, double, VectorizedArray<double>>(0);
test<2, 2, 3, double, VectorizedArray<double>>(1);
test<2, 2, 3, double, VectorizedArray<double>>(2);
}

0 comments on commit be8f0d4

Please sign in to comment.