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
seein 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<intdim>
classPointHistory
{ ...
// Fully referential configuration
SymmetricTensor<2, dim>
get_S() const
{
returnPhysics::Transformations::Contravariant::pull_back(get_tau(),
F);
}
SymmetricTensor<4, dim>
get_H() const
{
returnPhysics::Transformations::Contravariant::pull_back(get_Jc(),
F);
}
// Mixed / two-point configuration
Tensor<2, dim>
get_P() const
{
returnget_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
constSymmetricTensor<2, dim>S =get_S();
constSymmetricTensor<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(unsignedintI =0; I <dim; ++I)
for(unsignedintJ =0; J <dim; ++J)
for(unsignedintk =0; k <dim; ++k)
for(unsignedintL =0; L <dim; ++L)
for(unsignedintK =0; K <dim; ++K)
tmp1[I][J][k][L] +=get_F()[k][K] *H[I][J][K][L];
Tensor<4, dim>HH_mixed;
constTensor<2, dim>I_ns
=Physics::Elasticity::StandardTensors<dim>::I;
for(unsignedinti =0; i <dim; ++i)
for(unsignedintJ =0; J <dim; ++J)
for(unsignedintk =0; k <dim; ++k)
for(unsignedintL =0; L <dim; ++L)
{
// Push forward index 0 of material stress contribution
for(unsignedintI =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];
}
returnHH_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(unsignedintq_point =0; q_point <this->n_q_points;
++q_point) { const Tensor<2,dim> P = lqph[q_point]->get_P();
for(unsignedinti =0; i <this->dofs_per_cell; ++i) {
constunsignedinti_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(unsignedintq_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(unsignedinti =0; i <this->dofs_per_cell; ++i) {
constunsignedinti_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(unsignedintj =0; j <=i; ++j)
{
constunsignedintj_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
constTensor<2,dim,VectorizedArrayType>grad_u=
phi.get_gradient(q_point);
// Compute deformation gradient
constTensor<2,dim,VectorizedArrayType>F=
Physics::Elasticity::Kinematics::F(grad_u);
constSymmetricTensor<2,dim,VectorizedArrayType>b=
Physics::Elasticity::Kinematics::b(F);
constVectorizedArrayTypedet_F=determinant(F);
// Invert F
constTensor<2,dim,VectorizedArrayType>F_inv=invert(F);
constSymmetricTensor<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);
//Gradientitslefisincludedintheintegratecall, apply push
forward with F^{-1}
constTensor<2,dim,VectorizedArrayType>symm_grad_Nx=
symmetrize(F_inv);
// Compute the result
constTensor<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);
//endloopoverquadraturepoints
// 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/
<http://www.dealii.org/>
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
<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 [email protected].
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/
<http://www.dealii.org/>
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
<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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/1062eaa4-0e32-43f8-8fc4-70dfe2638a81n%40googlegroups.com
<https://groups.google.com/d/msgid/dealii/1062eaa4-0e32-43f8-8fc4-70dfe2638a81n%40googlegroups.com?utm_medium=email&utm_source=footer>.
