Hi Jean-Paul,
Thank you for your detailed response. I checked everything twice and found
my code consistent with what you suggested. I think I understood the issue
of Newton method convergence. The issue is occurring when I am giving a
very small compressive displacement as a boundary condition (near -1e-06).
Is there anyway I can get the solution converged even with a very small
compressive displacement?
Thanks!
Animesh
I am getting something like following -
Assembly method: Residual and linearisation are computed manually.
Grid:
Reference volume: 1200
Triangulation:
Number of active cells: 4800
Number ofset degrees of freedom: 9922
Setting up quadrature point data...
Timestep 1 @ 1s
_______________________________________________________________________________________
SOLVER STEP | LIN_IT LIN_RES RES_NORM RES_U
NU_NORM NU_U
_______________________________________________________________________________________
0 CST ASM SLV | 1 0.000e+00 1.000e+00 1.000e+00
1.000e+00 1.000e+00
1 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
2 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
3 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
4 CST ASM SLV | 1 0.000e+00 7.223e-09 7.223e-09
2.843e-08 2.843e-08
5 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
6 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
7 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
8 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
9 CST ASM SLV | 1 0.000e+00 7.223e-09 7.223e-09
2.843e-08 2.843e-08
10 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
11 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
12 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
13 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
14 CST ASM SLV | 1 0.000e+00 7.223e-09 7.223e-09
2.843e-08 2.843e-08
15 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
16 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
17 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
18 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
19 CST ASM SLV | 1 0.000e+00 7.223e-09 7.223e-09
2.843e-08 2.843e-08
20 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
21 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
22 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
23 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
24 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
25 CST ASM SLV | 1 0.000e+00 7.223e-09 7.223e-09
2.843e-08 2.843e-08
26 CST ASM SLV | 1 0.000e+00 1.696e-09 1.696e-09
6.675e-09 6.675e-09
..................
On Thursday, November 12, 2020 at 1:34:30 PM UTC+5:30 Jean-Paul Pelteret
wrote:
> Hi Animesh,
>
>
> > In the documentation in code gallery
> <https://dealii.org/developer/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html>,
>
> it is written that -
>
> This code is based off of step-44, and its documentation is largely
> identical to what's written there. It looks like that comment was missed
> when adjustments were being made to the documentation, so you should ignore
> it. You can see that the stress and tangent are computed directly using the
> kinematic variables passed as arguments to the class methods:
>
> const SymmetricTensor<2,dim,NumberType>
> <https://dealii.org/developer/doxygen/deal.II/classSymmetricTensor.html>
> tau = lqph[q_point]->get_tau(det_F,b_bar);
> const SymmetricTensor<4,dim,NumberType>
> <https://dealii.org/developer/doxygen/deal.II/classSymmetricTensor.html>
> Jc = lqph[q_point]->get_Jc(det_F,b_bar);
>
> > I would also need to update the quadrature point data to the previous
> step. Could you please let me know, how am I supposed to do that?
>
> You can see from the code snippet , as well as from the way that the
> assembly is structured
>
> template <int dim,typename NumberType>
> void Solid<dim,NumberType>::assemble_system(const BlockVector<double>
> <https://dealii.org/developer/doxygen/deal.II/classBlockVector.html>
> &solution_delta)
> {
> ...
> const BlockVector<double>
> <https://dealii.org/developer/doxygen/deal.II/classBlockVector.html>
> solution_total(get_total_solution(solution_delta));
> ...
> }
>
> template <int dim,typename NumberType>
> BlockVector<double>
> <https://dealii.org/developer/doxygen/deal.II/classBlockVector.html>
> Solid<dim,NumberType>::get_total_solution(const BlockVector<double>
> <https://dealii.org/developer/doxygen/deal.II/classBlockVector.html>
> &solution_delta) const
> {
> BlockVector<double>
> <https://dealii.org/developer/doxygen/deal.II/classBlockVector.html>
> solution_total(solution_n);
> solution_total += solution_delta;
> return solution_total;
> }
>
> that the only requirement to ensure a consistent state for the entire
> problem is that the addition of "solution_n" with "solution_delta" reflects
> the configuration of the body that you're wanting. So if that you're doing
> in the main "time loop" is something like
>
> while (time.current() <= time.end())
> {
> solution_delta = 0.0;
> solve_nonlinear_timestep(solution_delta);
> solution_n += solution_delta;
>
> auto tmp_soln_delta = solution_delta;
>
> solution_delta = 0.0; // get_total_solution() takes in the soltuion_delta,
> so don't solve the eigenvalue problem with the delta accounted for twice
> const bool accept_step = solve_eigenvalue_problem();
> if (accept_step == false)
> {
> solution_n -= tmp_soln_delta;
>
> adjust_boundary_conditions_or_timestep_size_to_apply_a_smaller_displacement_increment();
> }
>
> time.increment();
> }
>
> then that would seem a reasonable approach to me. As I said before, the
> stresses and strains are kept in sync with the solution_total within the
> assembly loop.
>
> I hope that this helps you to locate the source of your issues.
>
> Best,
> Jean-Paul
>
>
> On 10.11.20 16:06, Animesh Rastogi IIT Gandhinagar wrote:
>
> Hi Jean-Paul,
>
> Another question which might be the reason why I am getting the
> inconsistencies in my newton-raphson convergence. I just realised this..
>
> What I am doing is the following -
>
> I am giving small compressive displacements to the film+substrate. I am
> also calculating the eigenvalues of the tangent_stiffness matrix (after the
> convergence of newton method), to understand the stability of the system.
> The moment I encounter a negative smallest eigenvalue, I want to go to a
> previous time step and would like to solve for a small compressive
> displacement. For instance, if the compression is 0.1 at this time step and
> the eigenvalue is negative, I go to the previous time step and give a
> compressive displacement of 0.05.
>
> Here is my concern, I have added a condition that, if I find a negative
> eigenvalue during any time step, I update the *solution_n = solution_n -
> solution_delta*; and then run the time step again (with a smaller
> compressive displacement) after doing *solution_delta = 0.0*. Since the
> next state is also dependent on the updated quadrature point data (stresses
> etc.), I would also need to update the quadrature point data to the
> previous step. Could you please let me know, how am I supposed to do that?
>
> In the documentation in code gallery
> <https://dealii.org/developer/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html>,
>
> it is written that -
>
> "*The first function is used to create a material object and to
> initialize all tensors correctly: The second one updates the stored values
> and stresses based on the current deformation measure **Grad**u**n*"
> (inside the Point History class). However, I only find one function which
> is setup_lqp, and not the update_lqp or something similar.
>
> Thanks!
>
> Animesh
>
> On Tuesday, November 10, 2020 at 5:24:13 PM UTC+5:30 Animesh Rastogi IIT
> Gandhinagar wrote:
>
>> Hi Jean-Paul,
>>
>> There seem to be some inconsistencies in the simulations results that I
>> am getting, therefore I wanted to check if I have applied the material laws
>> properly accounting for the heterogeneities. Below, I have addted the
>> edited functions that are responsible for that. Could you please let me
>> know if there is something wrong here? I have checked the material
>> properties (mu1, mu2 etc), and they are being read correctly by the code
>> from the parameters file. I have also checked that the material ids (1 and
>> 2) are correctly applied to the cells where I need them.
>>
>> --------------------------------------
>> void setup_lqp (const Parameters::AllParameters ¶meters, const
>> unsigned int &materialid)
>> {
>>
>> //Material 1 - Thin Film
>> if (materialid == 1)
>> {
>>
>> std::cout<<"1_id"<<std::endl<<parameters.mu1<<std::endl<<parameters.nu1<<std::endl;
>> material.reset(new
>> Material_Compressible_Neo_Hook_One_Field<dim,NumberType>(parameters.mu1,
>> parameters.nu1));
>> }
>> //Material 2 - Substrate
>> if (materialid ==2)
>> {
>> material.reset(new
>> Material_Compressible_Neo_Hook_One_Field<dim,NumberType>(parameters.mu2,
>> parameters.nu2));
>> }
>>
>> }
>>
>> -------------------------------------------
>>
>> template <int dim,typename NumberType>
>> void Solid<dim,NumberType>::setup_qph()
>> {
>> std::cout << " Setting up quadrature point data..." << std::endl;
>>
>> quadrature_point_history.initialize(triangulation.begin_active(),
>> triangulation.end(),
>> n_q_points);
>>
>> for (typename Triangulation<dim>::active_cell_iterator cell =
>> triangulation.begin_active(); cell != triangulation.end();
>> ++cell)
>> {
>> const std::vector<std::shared_ptr<PointHistory<dim,NumberType> >
>> > lqph =
>> quadrature_point_history.get_data(cell);
>> Assert(lqph.size() == n_q_points, ExcInternalError());
>>
>> for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
>> lqph[q_point]->setup_lqp(parameters, cell->material_id());
>> }
>> }
>>
>> --------------------------------------
>>
>> Thanks!
>>
>> Animesh
>> On Tuesday, November 3, 2020 at 12:09:03 PM UTC+5:30 Animesh Rastogi IIT
>> Gandhinagar wrote:
>>
>>> Hi Jean,
>>>
>>> Sorry for late response. Thank you for your detailed answer. Using your
>>> suggestions I was able to assign the material ids to different area of my
>>> domain of interest.
>>>
>>> Thanks again!
>>>
>>> Animesh
>>>
>>> On Tuesday, October 27, 2020 at 3:05:39 AM UTC+5:30 Jean-Paul Pelteret
>>> wrote:
>>>
>>>> Dear Animesh,
>>>>
>>>> Let me preface my suggestion by saying that you can introduce some
>>>> heterogeneity in more than one way. That is, you're not bound to using a
>>>> material ID (which is what I'm going to suggest), but you could use some
>>>> geometric arguments to determine which class or material properties should
>>>> dictate the constitutive response at a continuum point. You also need not
>>>> use the class-based structure that this code galley example has. There are
>>>> other tutorials (e.g. step-8) that use Functions that return the material
>>>> coefficients at any point in space. So you should keep that in mind as you
>>>> read my suggestion.
>>>>
>>>> In my opinion, the simplest approach would be the following:
>>>>
>>>> 1. The PointHistory::setup_lqp() function
>>>> <https://github.com/dealii/code-gallery/blob/master/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc#L720-L724>
>>>>
>>>> is the point at which the constitutive law at some quadrature point is
>>>> set.
>>>> So you should extend this function in some way to make a decision about
>>>> which material it is in, and what the material law at that point is. The
>>>> first thing that you might do is add an argument for the material ID, and
>>>> then choose constiutitve parameters based on the material ID. You'll note
>>>> that the concrete constitutive law is stored in a pointer. This makes it
>>>> easy to extend to support other constitutive laws, because you can for
>>>> example abstract the Material_Compressible_Neo_Hook_One_Field class into
>>>> some base class and other classes (maybe a
>>>> Material_Compressible_Mooney_Rivlin_One_Field) and then, based on the
>>>> material ID, select which constitutive law is applied at that quadrature
>>>> point.
>>>> 2. The PointHistory::setup_lqp() function is called from
>>>> Solid::setup_qph()
>>>> <https://github.com/dealii/code-gallery/blob/master/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc#L1203-L1211>
>>>>
>>>> , and more specifically inside a loop over all cells. Its possible to
>>>> retrieve the cell->material_id() or cell->center() (or, if you create a
>>>> quadrature rule to help you then you can get the actual position of the
>>>> quadrature point if you want such granularity). You can then pass this
>>>> information in to PointHistory::setup_lqp() to help decide which
>>>> constitutive law / parameters to apply. If using the material ID then you
>>>> will, of course, want to set the material ID when building or reading in
>>>> your triangulation.
>>>>
>>>> That is, at least, how I typically do things when using this code as a
>>>> basis for extension. I hope that this helps you introduce heterogeneities
>>>> into your problem!
>>>>
>>>> Best,
>>>> Jean-Paul
>>>>
>>>>
>>>> On 26.10.20 08:15, Animesh Rastogi IIT Gandhinagar wrote:
>>>>
>>>> Hello All,
>>>>
>>>> I am trying to solve the wrinkling problem of a thin stiff film
>>>> attached to a soft substrate in dealii. I am using the code - from the
>>>> code
>>>> gallery for this purpose. For the simulations, I would have to consider
>>>> the
>>>> film and the substrate with different material properties.
>>>>
>>>> I understand that I would have to do something within the class
>>>> Qudrature point history. However, I am getting a little confused on how to
>>>> approach it. It would be great if someone could help me with this.
>>>>
>>>> This is the code that I am using for this purpose -
>>>> https://dealii.org/developer/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html
>>>>
>>>> Thanks!
>>>>
>>>> Animesh
>>>>
>>>> --
>>>> 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/26db8831-e190-4fcd-8247-b7172ab880fen%40googlegroups.com
>>>>
>>>> <https://groups.google.com/d/msgid/dealii/26db8831-e190-4fcd-8247-b7172ab880fen%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+un...@googlegroups.com.
>
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/f1209ac5-a4b8-46d4-84f3-12ca38f7d8f2n%40googlegroups.com
>
> <https://groups.google.com/d/msgid/dealii/f1209ac5-a4b8-46d4-84f3-12ca38f7d8f2n%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/4e377c7f-240f-4130-a82e-f2e77fee5fafn%40googlegroups.com.
