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

Matrix-free Piola transformation for affine cells #13642

Merged
merged 5 commits into from
May 23, 2022
+1,570 −813
Merged

Matrix-free Piola transformation for affine cells #13642

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5ac08a7
Piola transform for affine cells
NiklasWik May 2, 2022
ea402a8
Use ndarray instead of c-array
NiklasWik May 2, 2022
84a3be4
Avoid the use of determinant()
NiklasWik May 3, 2022
6f8c74d
Move template normal_dir from function to struct
NiklasWik May 9, 2022
41c9c46
Avoid copy paste
NiklasWik May 23, 2022
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
1,120 changes: 558 additions & 562 deletions include/deal.II/matrix_free/evaluation_kernels.h
View file
Open in desktop

Large diffs are not rendered by default.

489 changes: 444 additions & 45 deletions include/deal.II/matrix_free/fe_evaluation.h
View file
Open in desktop

Large diffs are not rendered by default.

12 changes: 10 additions & 2 deletions include/deal.II/matrix_free/fe_evaluation_data.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -735,8 +735,16 @@ 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 the first inverse transpose Jacobian (q_point = 0) is set for
* Cartesian/affine cells, and the actual Jacobian is stored at index 1
* instead. 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.
* Filled from MappingInfoStorage.jacobians in
* include/deal.II/matrix_free/mapping_info.templates.h
*/
const Tensor<2, dim, Number> *jacobian;

Expand Down
26 changes: 25 additions & 1 deletion include/deal.II/matrix_free/mapping_info.templates.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2177,6 +2177,18 @@ namespace internal
vv,
my_data.jacobians[0][offset + q][d][e]);
}
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);
Comment on lines +2183 to +2185
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You need here to pass in the reference-cell type. See other occurrences. Maybe we should remove the default parameter there so that such a mistake does not happen again?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I figured that I didn't have to since it is not done on line 2172. Is it a bug?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well it is. I wasn't sure which one of the two that was being used. Looking at it again makes it clear that I don't need this part, but only the same thing I did a couple of lines above. Since this is not used by the piola transform, and not anywhere else, I see no point in keeping it. I can keep it though, for consistency, and just remove the reorder_face_derivative_indices since it essentially does not do anything.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or I mean it does Assert(index < dim, ExcInternalError()); so I guess it has some value.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is actually only used in 2 out of 12 places with reference-cell passed in. And what I said above is not correct. Sorry

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just see this is the MappingQ path. Simplex meshes and mixed meshes should not enter this path; so you don't need to pass in the extra argument. Sorry for the confusion.

for (unsigned int d = 0; d < dim; ++d)
store_vectorized_array(
jac[d][ee],
vv,
my_data.jacobians[0][offset + q + 1][d][e]);
}

if (update_flags_faces & update_jacobian_grads)
{
Expand Down Expand Up @@ -2281,6 +2293,18 @@ namespace internal
vv,
my_data.jacobians[1][offset + q][d][e]);
}
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][d][e]);
}

if (update_flags_faces & update_jacobian_grads)
{
Expand Down Expand Up @@ -2762,7 +2786,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));
Comment on lines -2765 to +2789
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Was this a bug or is it related to the addition above?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The addition!

}

const UpdateFlags update_flags_common =
Expand Down
8 changes: 8 additions & 0 deletions include/deal.II/matrix_free/mapping_info_storage.h
View file
Open in desktop
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