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 dealii+unsubscr...@googlegroups.com
<mailto:dealii+unsubscr...@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/0c4f8f47-31ed-eddd-ac0d-d960c2c8c6cb%40gmail.com.
