-
Notifications
You must be signed in to change notification settings - Fork 708
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.
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
Matrix-free Piola transformation for affine cells #13642
Changes from 4 commits
5ac08a7
ea402a8
84a3be4
6f8c74d
41c9c46
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Original file line number | Diff line number | Diff line change | ||||
---|---|---|---|---|---|---|
|
@@ -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 set to a useful value if on a non-Cartesian cell. If the cell is | ||||||
* Cartesian/affine then the Jacobian is stored at index 1. For faces on | ||||||
* hypercube elements, the derivatives are reorder s.t the derivative | ||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Nit-pick:
Suggested change
|
||||||
* 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; | ||||||
|
||||||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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? There was a problem hiding this comment. Choose a reason for hiding this commentThe 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? There was a problem hiding this comment. Choose a reason for hiding this commentThe 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 There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Or I mean it does There was a problem hiding this comment. Choose a reason for hiding this commentThe 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 There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I just see this is the |
||
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) | ||
{ | ||
|
@@ -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) | ||
{ | ||
|
@@ -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
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Was this a bug or is it related to the addition above? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The addition! |
||
} | ||
|
||
const UpdateFlags update_flags_common = | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What do you mean by
Only set to a useful value if on a non-Cartesian cell.
? Do you mean that the index only goes by the quadrature pointsq
for the non-affine case, while for affine cells we always need to index with[0]
? And that the affine/Cartesian case also has the actual Jacobian in the same location?There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This was something that was already there and frankly something I too was confused about, but something I left untouched. I'll clearify it with that for affine/Cartesian cells only index 0 is set.