Matrix-free Piola transformation for affine cells #13642

NiklasWik · 2022-04-25T15:03:14Z

This PR builds on #13591 and adds the Piola transformation for affine cells when using the matrix-free evaluation of the FE_RaviartThomasNodal element. It also updates the tests associated with the previous pr (with ideas from @peterrum) and adds two new ones for non-Cartesian cells.

( #9545)

peterrum · 2022-04-25T20:51:51Z

@NiklasWik I'll take a look at this PR once the other PR is merged and you have rebased this one!

kronbichler · 2022-04-27T18:03:34Z

include/deal.II/matrix_free/fe_evaluation.h

+        {
+          const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
+            this->jacobian[0];
+          const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);


This is not documented anywhere, but you do not need to compute inverse matrix for Cartesian/affine cells

Suggested change

const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);

const Tensor<2, dim, VectorizedArrayType> &t_jac = this->jacobian[1];

This is filled in here for cells:

dealii/include/deal.II/matrix_free/mapping_info.templates.h

Lines 1212 to 1217 in 7e1892d

for (unsigned int d = 0; d < dim; ++d)

for (unsigned int e = 0; e < dim; ++e)

store_vectorized_array(

jac[d][e],

vv,

my_data.jacobians[0][idx + 1][d][e]);

. Note that this is not done for faces (if I remember correctly) but we could do it if necessary).

Ok! It does seem unnecessary to compute it twice... I'll have a look at it tomorrow (and also update the documentation)

I can't find that this is used anywhere in matrix_free/fe_evaluation.h but I'm unsure if it is used anywhere else. Because I would like to change it to the transpose jacobian instead of just the jacobian. That way there's no need to call transpose(invert(inv_t_jac)) if the cell is general in for instance get_values. I would also add it for faces.

The use is here:

dealii/include/deal.II/matrix_free/fe_evaluation_data.h

Lines 1276 to 1292 in c02086a

if (is_face == false &&

this->cell_type <= internal::MatrixFreeFunctions::affine)

{

Assert(this->jacobian != nullptr, ExcNotInitialized());

Point<dim, Number> point = this->quadrature_points[0];

const Tensor<2, dim, Number> &jac = this->jacobian[1];

if (this->cell_type == internal::MatrixFreeFunctions::cartesian)

for (unsigned int d = 0; d < dim; ++d)

point[d] += jac[d][d] * static_cast<typename Number::value_type>(

this->descriptor->quadrature.point(q)[d]);

else

for (unsigned int d = 0; d < dim; ++d)

for (unsigned int e = 0; e < dim; ++e)

point[d] += jac[d][e] * static_cast<typename Number::value_type>(

this->descriptor->quadrature.point(q)[e]);

return point;

Regarding a call to transpose, I would expect that an optimizing compiler optimizes the memory operations away. (It does not always work, but I've seen it work often.) Hence, it should not matter whether you call transpose(jac). Note that for general cells the Jacobians are not stored, just for the affine/Cartesian case (as the latter would be too memory intensive).

kronbichler · 2022-04-27T18:14:35Z

source/dofs/dof_tools_constraints.inst.in

+    template void DoFTools::make_periodicity_constraints(
+      const DoFHandler<deal_II_dimension,
+                       deal_II_space_dimension>::face_iterator &,
+      const DoFHandler<deal_II_dimension,
+                       deal_II_space_dimension>::face_iterator &,
+      AffineConstraints<float> &,
+      const ComponentMask &,
+      bool,
+      bool,
+      bool,
+      const FullMatrix<double> &,
+      const std::vector<unsigned int> &,
+      const float);
+


I think this change in the number of instantiations and the one a few lines below should be in a separate PR. Would you mind opening that one so it can be reviewed separately?

Yes, I agree!

NiklasWik · 2022-04-28T12:09:06Z

I ended up changing to the transpose jacobian anyway because it simplifies the code a bit. Also added it to the faces.

All matrixfree tests pass.

NiklasWik · 2022-04-28T14:01:32Z

I guess it is used somewhere else as well. I'll change it back, but the code will get a bit uglier. :/

kronbichler · 2022-04-29T09:26:58Z

/rebuild

peterrum

@NiklasWik Let me post what I have till now. I will continue with include/deal.II/matrix_free/fe_evaluation.h later today.

peterrum · 2022-04-29T09:14:56Z

tests/matrix_free/matrix_vector_rt_face_02.cc

+// This function tests the correctness of the matrix-free implementation
+// of the FE_RaviartThomasNodal element by evaluating values + gradients
+// as well as the divergence on faces and comparing the result with
+// FEFaceVaules which is considered the reference. The mesh is a
+// subdivided_parallelepiped with no hanging nodes and no other constraints,
+// i.e affine but not cartesian.


You could shorten the comment here by referring to the other test and only state what is different (i.e., non-Cartesian).

peterrum · 2022-04-29T09:15:33Z

tests/matrix_free/matrix_vector_rt_face_common.h

@@ -61,7 +61,38 @@ template <int dim, int fe_degree>
 void
 test();

+enum TestType : unsigned char


peterrum · 2022-04-29T09:16:55Z

tests/matrix_free/matrix_vector_rt_common.h

+std::string
+enum_to_string(TestType const enum_type)


A bit over-engineered :)

I just grabbed it from exadg I think

peterrum · 2022-04-29T09:18:55Z

include/deal.II/matrix_free/mapping_info.templates.h

-                              my_data.jacobians[0][idx + 1][d][e]);
+                              my_data.jacobians[0][idx + 1][e][d]);


Is this a bug? If yes, could you extract it to separate PR?

No it's not! It has to do with the discussion above with storing the transpose or not. Not sure if this is what breaks the tests or not. I made a change in include/deal.II/matrix_free/fe_evaluation_data.h where this is used.

peterrum · 2022-04-29T09:21:26Z

include/deal.II/matrix_free/mapping_info.templates.h

+                        const unsigned int ee =
+                          ExtractFaceHelper::reorder_face_derivative_indices<
+                            dim>(interior_face_no, e);


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?

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

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.

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

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

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.

peterrum · 2022-04-29T09:34:34Z

include/deal.II/matrix_free/fe_evaluation.h

+{
+  template <int dim, bool is_face>
+  inline const unsigned int *
+  get_index(const std::uint8_t fn)


What is fn?

Face number

peterrum · 2022-04-29T09:34:53Z

include/deal.II/matrix_free/fe_evaluation.h

+{
+  template <int dim, bool is_face>
+  inline const unsigned int *
+  get_index(const std::uint8_t fn)


What does this function do? Why does it return a pointer?

I use it to get the right indicies for faces, i.e the inverse of reorder_face_derivative_indices

peterrum · 2022-04-29T09:36:41Z

include/deal.II/matrix_free/fe_evaluation.h

+    static const unsigned int index0[3] = {0, 1, 2};
+    static const unsigned int index1[3] = {1, 0, 2};
+    static const unsigned int index2[3] = {2, 0, 1};
+    static const unsigned int index3[3] = {1, 2, 0};


Lets' try not use C-arrays and rather use std::array? I guess what you want here is 2D array? You can check out ndarray.

I didn't really want to call this function more than once, and I figured that I only wanted to return a 1D array. But maybe I should do like what is done in reorder_face_derivative_indices, and call it for every index?

I think that a static std::array should not be all that different from the pointer to a 1D array, as long as the compiler inlines everything correctly. I can't say whether it makes sense to call a function similar to reorder_face_derivative_indices for every index without benchmarking or looking at the assembly output of a compiler, as the function here is different from the one in MappingInfo in the sense that this one here is a hot loop, whereas the other function is "only" a setup function and thus potentially not used as often. That was also the intent for the original design by @kkormann and me, we tried to move everything expensive like shuffling around indices to the setup function to have as simple data structures in the actual operator evaluation for the regular DG face integrals.

Obviously, here we need to apply the reordering for every point in the evaluation. Or do we? Right now it is only curiosity, but do you think one could simplify the evaluators by shuffling around the normal and tangential components of the solution in such a way that the re-ordering of the Jacobian matrix on the faces would be the correct one again? I could even imagine that this might somewhat reduce the number of cases in the evaluation_kernels for the faces again, at least in terms of polynomial degrees that need to be expanded in normal and tangential direction. What do you think?

I do think you are correct! I see no immediate reason why we wouldn't be able to shuffle the components like we shuffle the derivatives in gradients_quad for instance, which is why we don't have to shuffle the inverse jacobian in the get/submit functions. I'll look into it.

Just fyi, it really made the function in evaluation_kernels muuuch nicer! Great idea @kronbichler, thanks! :)

Now I also (finally) see the value of having zx-coordinates for the faces

peterrum · 2022-04-29T09:39:07Z

include/deal.II/matrix_free/fe_evaluation.h

+  /**
+   * @copydoc FEEvaluationBase<dim,dim,Number,is_face>::get_value()
+   */
+  value_type
+  get_value(const unsigned int q_point) const;


Is this to make similar to what is done for get_gradient()?

peterrum · 2022-04-29T09:41:36Z

include/deal.II/matrix_free/fe_evaluation.h

+  if (this->data->element_type ==
+      internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
+    {
+      // TODO
+      AssertThrow(false, ExcNotImplemented());
+    }


Suggested change

if (this->data->element_type ==

internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)

{

// TODO

AssertThrow(false, ExcNotImplemented());

}

AssertThrow(this->data->element_type !=

internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas, ExcNotImplemented());

The if-branch will be needed later when this gets implemented anyway. Or are there any other benefits other than being more nicer looking atm?

peterrum

Looks good to me! If the simplex tests should pass, what is missing in this PR? The general case can be done in a follow-up PR!?

peterrum · 2022-05-02T09:28:13Z

include/deal.II/matrix_free/mapping_info.templates.h

+                        const unsigned int ee =
+                          ExtractFaceHelper::reorder_face_derivative_indices<
+                            dim>(interior_face_no, e);


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.

peterrum · 2022-05-02T09:36:03Z

include/deal.II/matrix_free/evaluation_kernels.h

+      if (dim == 3)
+        evaluate_in_face_apply<1, EvalTangent, EvalNormal>(


Why is the if-statement for dim==3 and not below for the direction 2 case?

The order here shouldn't actually matter because I do all the offsets in evaluate_in_face_apply. However this makes more sense since I now always put the evaluation of the normal faces last, i.e face_direction == component . It also made it easier for me to write everything up on paper beforehand.

But what is face_direction==2 in the 2D case?

Ah, the different calls are not with different face_direction but rather the normal_dir that I pass to the tensor_evaluation_kernels. Since I'm using the "general" kernels for the normal faces, this is not used for anything other than the initial offsets of the pointers.

Let me reformulate the question: but what is normal_dir==2 in the 2D case ;)

What about swapping these?

const unsigned int face_direction = face_no / 2; evaluate_in_face_apply<0, EvalGeneral, EvalGeneral>( values_dofs, fe_eval, scratch_data, evaluation_flag, face_direction, subface_index, std::integral_constant<bool, integrate>()); // normal evaluate_in_face_apply<1, EvalNormal, EvalTangent>( values_dofs, fe_eval, scratch_data, evaluation_flag, face_direction, subface_index, std::integral_constant<bool, integrate>()); // tangential 0 if (dim == 3) evaluate_in_face_apply<2, EvalTangent, EvalNormal>( values_dofs, fe_eval, scratch_data, evaluation_flag, face_direction, subface_index, std::integral_constant<bool, integrate>()); // tangential 1

Here

dealii/include/deal.II/matrix_free/evaluation_kernels.h

Lines 3190 to 3198 in 0557c00

case 3:

eval0.template values<0, true, false, normal_dir>(values_dofs,

values_quad);

eval1.template values<1, true, false, normal_dir>(values_quad,

values_quad);

break;

case 2:

eval0.template values<0, true, false, normal_dir>(values_dofs,

values_quad);

normal_dir can be 2 although the dimension of EvaluatorTensorProduct is only 1 or 2...

I see you have introduced normal_dir in EvaluatorTensorProduct<evaluate_general> during the commit NiklasWik@6918f72#diff-bd08275179f5b90beff71f4688219f973ca2cb77a58e1df1cb809aff9274fce1. I have overlooked that. Could you revert that change? There is no normal-direction in the general case. I guess the only reason for that was to be able to pass normal_dir to the RT type of EvaluatorTensorProduct. This does not feel right, since we might need in the future also something else than the normal. I would suggest to encode the information of the normal in the template arguments of the class.

What about swapping these?

I prefer this order since this is the order of the faces in the quadrature points. I think that changing it only will introduce more confusion. Unless it's beneficial from a performance pov, which I doubt..? I can add a comment explaining the order if you want.

Could you revert that change?

This will require a lot of changes... And I'm not immediately sure how I would make this pretty. I would rather skip this for now, and rather address this if I have time over at the end of my thesis. But if you really think this is necessary now then I'll of course do it!

normal_dir can be 2 although the dimension of EvaluatorTensorProduct is only 1 or 2...

Maybe my explanation above wasn't that good :')

peterrum · 2022-05-02T09:39:55Z

include/deal.II/matrix_free/evaluation_kernels.h

+      const unsigned int component_table[3][3] = {{1, 2, 0},
+                                                  {2, 0, 1},
+                                                  {0, 1, 2}};


Could you use ndarray? Also I think this could be static.

peterrum · 2022-05-02T09:41:36Z

include/deal.II/matrix_free/fe_evaluation.h

+          for (unsigned int comp = 0; comp < n_components; ++comp)
+            value_out[comp] = this->values_quad[comp * nqp + q_point] *
+                              jac[comp][comp] *
+                              inv_det; // / this->jacobian[0][comp][comp];


Is the comment still needed?

Not really. I was curious about whether or not it would be faster to use division here than to compute the determinant. @kronbichler thought it would be faster with the determinant (and I of course trust him on this) but I haven't measured it.

peterrum · 2022-05-02T09:44:04Z

include/deal.II/matrix_free/fe_evaluation.h

+          // Derivatives are reordered for faces. Need to take this into account
+          const VectorizedArrayType inv_det =
+            (is_face && dim == 2 && this->get_face_no() < 2) ?
+              -determinant(inv_t_jac) :
+              determinant(inv_t_jac);


Nope! The determinant is the same for all cases in 3D :)

peterrum · 2022-05-02T09:51:01Z

include/deal.II/matrix_free/fe_evaluation.h

+        {
+          const Tensor<2, dim, dealii::VectorizedArray<Number>> jac =
+            this->jacobian[1];
+          const VectorizedArrayType inv_det = determinant(this->jacobian[0]);


Shouldn't this be simply the product of the values on the diagonal?

Oh yes, sorry for the confusion. The division was instead of previously computing the inverse of this->jacobian[0]. But yes, maybe the product of the diagonal elements is faster than computing the determinant. What do you think?

Yes indeed, the product of the diagonals is faster (2 multiplications in 3D) compared to the determinant, which does 4 multiplications, 5 FMAs in 3D according to

dealii/include/deal.II/base/tensor.h

Lines 2749 to 2756 in 82cf0e0

// hard-coded for efficiency reasons

const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -

internal::NumberType<Number>::value(t[1][2] * t[2][1]);

const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -

internal::NumberType<Number>::value(t[1][0] * t[2][2]);

const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -

internal::NumberType<Number>::value(t[1][1] * t[2][0]);

return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;

peterrum · 2022-05-02T09:51:57Z

include/deal.II/matrix_free/fe_evaluation.h

+              determinant(inv_t_jac);
+
+          VectorizedArrayType tmp;
+          // J * grad_quad * J^-1 * det(J^-1)


Could you also add the equation to the Cartesian case?

peterrum · 2022-05-02T09:58:38Z

include/deal.II/matrix_free/fe_evaluation.h

@@ -5945,6 +6373,12 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
    const SymmetricTensor<2, dim, VectorizedArrayType> sym_grad,
    const unsigned int                                 q_point)
 {
+  // TODO


Suggested change

// TODO

peterrum · 2022-05-02T10:00:21Z

include/deal.II/matrix_free/fe_evaluation.h

+      const std::size_t                   nqp = this->n_quadrature_points;
+      Tensor<1, dim, VectorizedArrayType> value_out;
+
+      // Cartesian cell


Personally, I would prefer if you'd move the comments inside the if-statement body.

peterrum · 2022-05-02T10:06:55Z

include/deal.II/matrix_free/fe_evaluation_data.h

+   * 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 transpose Jacobian is stored at index 1. For
+   * faces, the derivatives are reorder s.t the normal derivative is stored
+   * last, i.e for 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].


Does the documentation of MappingInfo have this information? If no could you add it there as well. Could you reference here the class/field. The documentation might get out of sync at some time.

peterrum

Good to go from my side once you have addressed https://github.com/dealii/dealii/pull/13642/files#r862829918 and solved my confusion!

NiklasWik · 2022-05-04T15:02:18Z

I'm pretty sure that it is the matrix assembly that is taking so long for the tests (on the topic of why we can't run the testcase with dim = degree = 3). Noticed this as I'm comparing the performance, and where the actual matrix-vector product takes a fraction of the total runtime. Still of course slow compared to the mf version.

NiklasWik · 2022-05-09T12:00:29Z

I moved the template parameter, but there are some changes in evaluation_kernels.h that I don't like, for instance in interpolate_generic_raviart_thomas_apply_face(..). There's probably a better way to do this.

Also added a little comment with a table regarding the ordering of the faces, and made the small changes from #13671 for hanging nodes support. There's ofc still the problem with the tests...

Updated tests, reordering of components in face evaluation, and storing jacobian for affine face evaluation

Update comments and documentation

+ Bug fix for general submit_values

peterrum

Thanks for introducing EvaluatorTensorProductAnisotropic! I like it like this more! I guess by now you know me good enough that I would have preferred this change in a separate PR, since this change modifies a lot of unrelated places ;)

@kronbichler Could you take another look at the PR and after that let's merge!

include/deal.II/matrix_free/evaluation_kernels.h

peterrum · 2022-05-12T21:38:15Z

include/deal.II/matrix_free/evaluation_kernels.h

+      // -----------------------------------------------------------------------------------
+      // |          |                   Anisotropic faces                 | Isotropic faces|
+      // | Face dir | comp, coords, normal_dir | comp, coords, normal_dir | comp, coords   |
+      // | --------------------------------------------------------------------------------|


Thanks for adding a table! I am still not 100% happy with the code: it does not match the table. 2D does not contain any 2s; you remove the 2 division. I would suggest to get it on master, as it is and we clean up it later on.

I left out normal_dir for the isotropic faces in the table since the parameter is not used as normal_dir for those cases, but only for shifting the pointers. I'm not sure how you would like it to be instead.

Creating new EvaluationTensorProductAnisotropic

kronbichler · 2022-05-20T18:50:49Z

include/deal.II/matrix_free/fe_evaluation.h

+  if (this->data->element_type ==
+      internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
+    {
+      // Piola transform is required
+
+#  ifdef DEBUG
+      Assert(this->is_reinitialized, ExcNotInitialized());
+#  endif
+      AssertIndexRange(q_point, this->n_quadrature_points);
+      Assert(this->J_value != nullptr,
+             internal::ExcMatrixFreeAccessToUninitializedMappingField(
+               "update_gradients"));
+      Assert(this->jacobian != nullptr,
+             internal::ExcMatrixFreeAccessToUninitializedMappingField(
+               "update_gradients"));
+#  ifdef DEBUG
+      this->gradients_quad_submitted = true;
+#  endif
+
+      const std::size_t nqp = this->n_quadrature_points;
+      if (!is_face &&
+          this->cell_type == internal::MatrixFreeFunctions::cartesian)
+        {
+          // Cartesian cell
+          const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
+            this->jacobian[0];
+          const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
+          const VectorizedArrayType weight = this->quadrature_weights[q_point];
+          for (unsigned int d = 0; d < dim; ++d)
+            for (unsigned int comp = 0; comp < n_components; ++comp)
+              this->gradients_quad[(comp * dim + d) * nqp + q_point] =


Instead of copy-pasting this code, can we simply convert the Tensor<1, dim, Tensor<1, dim, VA>> into a Tensor<2, dim, VA>? Not sure if there is an automatic conversion, but otherwise you can simply copy the tensor and insert the code. It would be pity to repeat all the code here.

kronbichler · 2022-05-20T18:52:57Z

include/deal.II/matrix_free/fe_evaluation_data.h

-   * 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


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 points q 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?

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.

kronbichler · 2022-05-20T18:54:13Z

include/deal.II/matrix_free/fe_evaluation_data.h

+   * 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


Nit-pick:

Suggested change

* hypercube elements, the derivatives are reorder s.t the derivative

* hypercube elements, the derivatives are reorder s.t. the derivative

kronbichler

This is great work, it is clear that you put a lot of patience and effort into this work. I like the direction of this very much, we are almost there for the RT case.

Fix documentation

Matrix-free Piola transformation for affine cells

peterrum self-requested a review April 25, 2022 20:51

NiklasWik force-pushed the piola_matrix_free branch 2 times, most recently from 7f770f4 to 4632283 Compare April 27, 2022 08:35

kronbichler reviewed Apr 27, 2022

View reviewed changes

NiklasWik force-pushed the piola_matrix_free branch from 4632283 to 2dacda6 Compare April 28, 2022 12:06

kronbichler added Matrix-free ready to test labels Apr 29, 2022

peterrum reviewed Apr 29, 2022

View reviewed changes

NiklasWik force-pushed the piola_matrix_free branch from c8a1154 to 0550f16 Compare May 2, 2022 08:56

peterrum reviewed May 2, 2022

View reviewed changes

NiklasWik force-pushed the piola_matrix_free branch from 0550f16 to a3b2bce Compare May 2, 2022 14:29

peterrum approved these changes May 2, 2022

View reviewed changes

NiklasWik force-pushed the piola_matrix_free branch from 10091e5 to 14881f0 Compare May 9, 2022 11:01

NiklasWik added 3 commits May 12, 2022 10:58

Piola transform for affine cells

5ac08a7

Updated tests, reordering of components in face evaluation, and storing jacobian for affine face evaluation

Use ndarray instead of c-array

ea402a8

Update comments and documentation

Avoid the use of determinant()

84a3be4

+ Bug fix for general submit_values

NiklasWik force-pushed the piola_matrix_free branch 2 times, most recently from 750c2f2 to c529a0e Compare May 12, 2022 12:18

peterrum approved these changes May 12, 2022

View reviewed changes

Move template normal_dir from function to struct

6f8c74d

Creating new EvaluationTensorProductAnisotropic

NiklasWik force-pushed the piola_matrix_free branch from 2b7e176 to 6f8c74d Compare May 13, 2022 13:01

kronbichler reviewed May 20, 2022

View reviewed changes

kronbichler approved these changes May 20, 2022

View reviewed changes

Avoid copy paste

41c9c46

Fix documentation

kronbichler approved these changes May 23, 2022

View reviewed changes

kronbichler merged commit fac6d58 into dealii:master May 23, 2022

This was referenced Jun 16, 2022

Tests for hanging nodes mf RaviartThomas #13989

Merged

Non-affine piola transformation for mf RaviartThomas #13990

Merged

mkghadban pushed a commit to OpenFCST/dealii that referenced this pull request Sep 8, 2022

Merge pull request dealii#13642 from NiklasWik/piola_matrix_free

b389e3f

Matrix-free Piola transformation for affine cells

	const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
	const Tensor<2, dim, VectorizedArrayType> &t_jac = this->jacobian[1];

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

	if (is_face == false &&
	this->cell_type <= internal::MatrixFreeFunctions::affine)
	{
	Assert(this->jacobian != nullptr, ExcNotInitialized());
	Point<dim, Number> point = this->quadrature_points[0];

	const Tensor<2, dim, Number> &jac = this->jacobian[1];
	if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
	for (unsigned int d = 0; d < dim; ++d)
	point[d] += jac[d][d] * static_cast<typename Number::value_type>(
	this->descriptor->quadrature.point(q)[d]);
	else
	for (unsigned int d = 0; d < dim; ++d)
	for (unsigned int e = 0; e < dim; ++e)
	point[d] += jac[d][e] * static_cast<typename Number::value_type>(
	this->descriptor->quadrature.point(q)[e]);
	return point;

		my_data.jacobians[0][idx + 1][d][e]);
		my_data.jacobians[0][idx + 1][e][d]);

		if (dim == 3)
		evaluate_in_face_apply<1, EvalTangent, EvalNormal>(

	case 3:
	eval0.template values<0, true, false, normal_dir>(values_dofs,
	values_quad);
	eval1.template values<1, true, false, normal_dir>(values_quad,
	values_quad);
	break;
	case 2:
	eval0.template values<0, true, false, normal_dir>(values_dofs,
	values_quad);

	// hard-coded for efficiency reasons
	const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
	internal::NumberType<Number>::value(t[1][2] * t[2][1]);
	const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
	internal::NumberType<Number>::value(t[1][0] * t[2][2]);
	const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
	internal::NumberType<Number>::value(t[1][1] * t[2][0]);
	return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;

	* hypercube elements, the derivatives are reorder s.t the derivative
	* hypercube elements, the derivatives are reorder s.t. the derivative

		std::string
		enum_to_string(TestType const enum_type)

Matrix-free Piola transformation for affine cells #13642

Matrix-free Piola transformation for affine cells #13642

Conversation

NiklasWik commented Apr 25, 2022 • edited

peterrum commented Apr 25, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

NiklasWik Apr 28, 2022 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

NiklasWik commented Apr 28, 2022

NiklasWik commented Apr 28, 2022 • edited

kronbichler commented Apr 29, 2022

peterrum left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum May 4, 2022 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum left a comment

Choose a reason for hiding this comment

NiklasWik commented May 4, 2022

NiklasWik commented May 9, 2022 • edited

peterrum left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kronbichler left a comment

NiklasWik commented Apr 25, 2022 •

edited

NiklasWik Apr 28, 2022 •

edited

NiklasWik commented Apr 28, 2022 •

edited

peterrum May 4, 2022 •

edited

NiklasWik commented May 9, 2022 •

edited