Hi Jean-Paul, first of all, thank you very much for your comprehensive answer. It's a little bit unfortunate that there is no way for the deformed test gradient evaluation using matrix-free. However, your code snippets and the explanation are traceable and the required changes for the tangent operator seem acceptable (assuming the indexing in the get_HH() function does the right thing, otherwise it's probably not acceptable). As always, I was a little bit too ambitious with my time schedule for this task. I will document your suggestions in an issue and try to implement the rephrased problem in the medium run. In the end, the change is entirely performance motivated, since the residual assembly (without matrx-free) is comparatively expensive considering that it is only performed once per nonlinear iteration.
As an interesting side note: If I remove my wrong 'symm_grad_Nx' contribution from my shared code snippet, the remaining code corresponds already to a rephrased formulation using the Kirchoff stress and deformation gradient (as you have already noticed). Using this residual in combination with the "tau-Jc" tangent from step-44 (one-field adopted) leads to a convergence of the nonlinear scheme. So, it seems the linearization is still "correct enough" for the modified (two-point) frame of the residual. I will -of course- not rely on this half baked solution. Best regards, David Jean-Paul Pelteret schrieb am Dienstag, 29. Dezember 2020 um 20:06:17 UTC+1: > Hi David, > > So as I'm sure that you know, if you want to assemble the linear system > for quasi-static non-linear (hyper)elasticity with a referential DoFHandler > (using no Eulerian mapping to transform the shape functions to the spatial > setting) AND without doing a similar transformation by hand (as step-44 > does) then you are pretty much limited to using one of two > parameterisations. Either you have to choose the (two-point) Kirchoff > stress (P) and deformation gradient (F) as the work conjugate pairs, or the > (fully-referential) Piola-Kirchhoff stress (S) and the Green-Lagrange > strain tensor (E = 0.5(C-I) with C being the right Cauchy-Green tensor). > > Thinking about just the residual term: to the best of my knowledge (and > from what I recall from the work that we did here > <https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6336>), the > variation dC = symm(F^{T}.dF) = symm(F^{T} d(Grad(u))) -- with u > representing the displacement -- cannot be framed a way that the > matrix-free framework supports: the test function cannot be pre-multiplied > by a field variable. One would have to somehow rephrase the whole problem > such that you end up with just the test function d(Grad(u)) as a > pre-multiplication to your stress variable, and in the end you come up with > the exactly the other parameterisation because P = F.S and dF = d(I + > Grad(u)) == d(Grad(u)). So you may as well start off by parameterising the > problem in terms of F and P, and naturally you have to adjust the > linearisation as well. The linearisation of the two-point tensor P contains > both the material and geometric tangent (you can identify each by > linearising the decomposition P = F.S to get these two quantities). I get > the sense from what I see in the code that you've shared you're trying to > accomplish exactly this (although I'm not sure its quite right, since tau = > P.F^{T} --> P = tau.F^{-T}, and you seem to have some extra manipulation > involving what you call symm_grad_Nx -- this looks a bit suspect to me). > > I have some code to share that will do the transformation from the stress > and tangent quantities already used in step-44 to those that you need for > this F-P parameterisation. You can put these in the PointHistory class > > template <int dim> > class PointHistory > { ... > // Fully referential configuration > SymmetricTensor<2, dim> > get_S() const > { > return Physics::Transformations::Contravariant::pull_back(get_tau(), F); > } > SymmetricTensor<4, dim> > get_H() const > { > return Physics::Transformations::Contravariant::pull_back(get_Jc(), F); > } > // Mixed / two-point configuration > Tensor<2, dim> > get_P() const > { > return get_F() * Tensor<2, dim>(get_S()); > } > Tensor<4, dim> > get_HH() const > { > // Reference: Simo1984a eq 4.4 // Simo, J. C. and Marsden, J. // On the > rotated stress tensor and the material version of the {Doyle-Ericksen} > formula // DOI: 10.1007/BF00281556 > const SymmetricTensor<2, dim> S = get_S(); > const SymmetricTensor<4, dim> H = get_H(); > // Push forward index 2 of material stress contribution > // This index operation relies on the symmetry of the > // last two indices, so therefore (in case it makes a > // difference) we transform it first. > Tensor<4, dim> tmp1; > for (unsigned int I = 0; I < dim; ++I) > for (unsigned int J = 0; J < dim; ++J) > for (unsigned int k = 0; k < dim; ++k) > for (unsigned int L = 0; L < dim; ++L) > for (unsigned int K = 0; K < dim; ++K) > tmp1[I][J][k][L] += get_F()[k][K] * H[I][J][K][L]; > Tensor<4, dim> HH_mixed; > const Tensor<2, dim> I_ns = Physics::Elasticity::StandardTensors<dim>::I; > for (unsigned int i = 0; i < dim; ++i) > for (unsigned int J = 0; J < dim; ++J) > for (unsigned int k = 0; k < dim; ++k) > for (unsigned int L = 0; L < dim; ++L) > { > // Push forward index 0 of material stress contribution > for (unsigned int I = 0; I < dim; ++I) > HH_mixed[i][J][k][L] += get_F()[i][I] * tmp1[I][J][k][L]; > // Add geometric stress contribution > HH_mixed[i][J][k][L] += I_ns[i][k] * S[J][L]; > } > return HH_mixed; > } > } > Naturally, the above could be simplified a bit for this fixed > parameterisation. > > With this, the (matrix-based) assembly would looks something like this for > the RHS > // Excuse the messiness -- coding in an email client is not a good idea! > for (unsigned int q_point = 0; q_point < this->n_q_points; ++q_point) { const > Tensor<2,dim> P = lqph[q_point]->get_P(); > for (unsigned int i = 0; i < this->dofs_per_cell; ++i) { > const unsigned int i_group = this->fe.system_to_base_index(i).first.first; > const Tensor<2,dim> dF_I = fe_values_ref[this->u_fe].gradient(i, q_point); > > if (i_group == this->u_dof) > data.cell_rhs(i) -= double_contract(dF_I, P) * JxW; > ... // See the rest of step-44 > > and the salient part of the tangent matrix assembly would be something like > for (unsigned int q_point = 0; q_point < this->n_q_points; ++q_point) { // > Linearisation of P with respect to F; // contains both material and > geometric tangent contributions const Tensor<4,dim> HH = > lqph[q_point]->get_HH(); > > for (unsigned int i = 0; i < this->dofs_per_cell; ++i) { > const unsigned int i_group = this->fe.system_to_base_index(i).first.first; > const Tensor<2,dim> dF_I = fe_values_ref[this->u_fe].gradient(i, q_point); > for (unsigned int j = 0; j <= i; ++j) > { > const unsigned int j_group = this->fe.system_to_base_index(j).first.first; > const Tensor<2,dim> dF_J = fe_values_ref[this->u_fe].gradient(j, q_point); > > if ((i_group == j_group) && (i_group == this->u_dof)) > { > cell_matrix(i, j) += scalar_product(dF_I, HH, dF_J) * JxW; } ... // See > the rest of step-44 > > So IIRC the call to phi.submit_gradient() that is bound to the > undeformed/non-mapped DoFHandler is effectively the same as testing against > what I've called dF_I in the above code. I think that the transformations > from tau->P and Jc->HH should make it simple enough to adjust what remains. > > I hope that this helps you implement what you're trying to do. > Jean-Paul > > > > On 29.12.20 15:59, 'David' via deal.II User Group wrote: > > Maybe as an edit: what I currently do looks the following way: > > ``` > // Get gradient in reference frame > const Tensor<2, dim, VectorizedArrayType> grad_u = phi.get_gradient( > q_point); > // Compute deformation gradient > const Tensor<2, dim, VectorizedArrayType> F = Physics::Elasticity:: > Kinematics::F(grad_u); > const SymmetricTensor<2, dim, VectorizedArrayType> b = Physics::Elasticity > ::Kinematics::b(F); > > const VectorizedArrayType det_F = determinant(F); > // Invert F > const Tensor<2, dim, VectorizedArrayType> F_inv = invert(F); > const SymmetricTensor<2, dim, VectorizedArrayType> b_bar = Physics:: > Elasticity::Kinematics::b( Physics::Elasticity::Kinematics::F_iso(F)); > // Get tau from material description > SymmetricTensor<2, dim, VectorizedArrayType> tau; material->get_tau(tau, > det_F, b_bar, b); > > // Gradient itslef is included in the integrate call, apply push forward > with F^{-1} > const Tensor<2, dim, VectorizedArrayType> symm_grad_Nx = symmetrize(F_inv > ); > // Compute the result > const Tensor<2, dim, VectorizedArrayType> result = symm_grad_Nx * Tensor<2 > , dim, VectorizedArrayType>(tau); > // Apply an additional push forward with F^{-T} > phi.submit_gradient(-result * transpose(F_inv), q_point); // end loop > over quadrature points > > // For each element > phi.integrate(false, true); > ``` > > It works, but I don't get the quadratic convergence order of the Newton > scheme anymore. Hence, I assume this is not sufficient. > David schrieb am Dienstag, 29. Dezember 2020 um 12:41:35 UTC+1: > >> >> Hi all, >> >> I'm currently trying to implement a vectorized variant of the residual >> assembly as it is done in step-44/one of the corresponding code-gallery >> examples using FEEvaluation objects in combination with matrix-free. >> >> I was not able to find an appropriate solution for the given code line >> <https://github.com/dealii/dealii/blob/57ca3f894d1bb31f6a6dd448864f4386dc6692ad/examples/step-44/step-44.cc#L2040>and >> >> the subsequent application for the assembly process: the major problem is >> that the shape function gradients are defined with regard to the current >> configuration, whereas my matrix-free object holds a mapping, which refers >> to the initial configuration. The reference configuration mapping of my >> matrix-free object is desired as the final integration is actually >> performed in the reference frame. >> A valid option is probably (I have not tried it) to use a matrix-free >> object with a different mapping (MappingQEulerian) which directly evaluates >> the shape gradients in the deformed configuration and use another >> referential matrix-free object for the integration part. >> The main reason to avoid this approach is that it would require a >> reinitialization of the 'deformed' matrix-free object in each nonlinear >> step and the mapping needs to be recomputed in each reinizialization. On >> the other hand, the expensive reinitialization is not really required as >> the mapping between the referential and the current configuration is known >> due to the (inverse) deformation gradient. >> Therefore, I'm looking for an opportunity to access the shape gradients >> and perform a push forward similarly to the approach in step-44 in order to >> evaluate the desired quantities. >> >> Until now I was not able to find an obvious solution for this >> computations using the FEEvaluation class. Maybe anyone else has a better >> idea for the described problem. >> >> Thanks in advance, >> David >> > -- > The deal.II project is located at http://www.dealii.org/ > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > --- > You received this message because you are subscribed to the Google Groups > "deal.II User Group" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to dealii+un...@googlegroups.com. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/7b3c1488-cf75-4526-8442-d1045dc73398n%40googlegroups.com > > <https://groups.google.com/d/msgid/dealii/7b3c1488-cf75-4526-8442-d1045dc73398n%40googlegroups.com?utm_medium=email&utm_source=footer> > . > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/1062eaa4-0e32-43f8-8fc4-70dfe2638a81n%40googlegroups.com.
