running the "The 'Goal-oriented mesh adaptivity in elastoplasticity problems' code gallery program "

[Muhammad Mashhood]

Dear deal.ii users,

                            I recently started working with deal.ii and want to work with material models from material mechanics. For this so far I came across nicely written and well worked code of "The 'Goal-oriented mesh adaptivity in elastoplasticity problems " in the code gallery. I tried to compile and run it on deal.ii 9.1.1 and one of the case ran but not giving displacement and stresses values in the loading step output results.

Anyone worked with it is warmly welcome to suggest or share any possible problem which I might not be taking care of? Thank you in advance!

Regards,

Muhammad

[Wolfgang Bangerth]

On 6/17/19 1:50 AM, Muhammad Mashhood wrote:
>                             I recently started working with deal.ii and want
> to work with material models from material mechanics. For this so far I came
> across nicely written and well worked code of "The 'Goal-oriented mesh
> adaptivity in elastoplasticity problems " in the code gallery. I tried to
> compile and run it on deal.ii 9.1.1 and one of the case ran but not giving
> displacement and stresses values in the loading step output results.

Can you be more specific what doesn't work? Which case is it that doesn't
work? And what do you mean by "not giving displacement and stresses in the ...
output"?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Muhammad Mashhood]

Thank you Prof. Wolfgang for the response! I will try my best to elaborate it more.

I ran the case of "Cantiliver_beam_3d" and it returned the results after couple of minutes as mentioned in the attachement (where per time step, seemingly only one iteration was performed in Newton method and in all increamental load steps the result values are zero for both stress and displacement).

Secondly I tried to run the case "Thick_tube_internal_pressure" which first popped up the dialog message of :

TimerOutput objects finalize timed values printed to the
screen by communicating over MPI in their destructors.
Since an exception is currently uncaught, this
synchronization (and subsequent output) will be skipped to
avoid a possible deadlock.
---------------------------------------------------------
ERROR: Uncaught exception in MPI_InitFinalize on proc 0. Skipping MPI_Finalize() to avoid a deadlock.


----------------------------------------------------
Exception on processing:

--------------------------------------------------------
An error occurred in line <3617> of file </home/muhammad/installed_software/dealii-9.1.1/examples/Goal-oriented_mesh_adaptivity_in_elastoplasticity_problems/elastoplastic.cc> in function
    void ElastoPlastic::ElastoPlasticProblem<dim>::make_grid() [with int dim = 3]
The violated condition was:
    dim == 2
Additional information:
    You are trying to use functionality in deal.II that is currently not implemented. In many cases, this indicates that there simply didn't appear much of a need for it, or that the author of the original code did not have the time to implement a particular case. If you hit this exception, it is therefore worth the time to look into the code to find out whether you may be able to implement the missing functionality. If you do, please consider providing a patch to the deal.II development sources (see the deal.II website on how to contribute).
--------------------------------------------------------

Aborting!

Where as per my understanding of template based nature of deall.ii and also this case is 2d so if I change the dim = 2 from originally dim = 3 in main function of elastoplastic.cc then the following error comes during compilation :

/home/muhammad/installed_software/dealii-9.1.1/examples/Goal-oriented_mesh_adaptivity_in_elastoplasticity_problems/elastoplastic.cc: In instantiation of ‘void ElastoPlastic::ElastoPlasticProblem<dim>::make_grid() [with int dim = 2]’:
/home/muhammad/installed_software/dealii-9.1.1/examples/Goal-oriented_mesh_adaptivity_in_elastoplasticity_problems/elastoplastic.cc:7100:14:   required from ‘void ElastoPlastic::ElastoPlasticProblem<dim>::run() [with int dim = 2]’
/home/muhammad/installed_software/dealii-9.1.1/examples/Goal-oriented_mesh_adaptivity_in_elastoplasticity_problems/elastoplastic.cc:7265:21:   required from here
/home/muhammad/installed_software/dealii-9.1.1/examples/Goal-oriented_mesh_adaptivity_in_elastoplasticity_problems/elastoplastic.cc:3796:30: error: invalid initialization of reference of type ‘dealii::Triangulation<3, 3>&’ from expression of type ‘dealii::parallel::distributed::Triangulation<2, 2>’
         extrude_triangulation(triangulation_2d,
                              ^
/home/muhammad/installed_software/dealii-9.1.1/examples/Goal-oriented_mesh_adaptivity_in_elastoplasticity_problems/elastoplastic.cc:107:3: note: in passing argument 4 of ‘void ElastoPlastic::extrude_triangulation(const dealii::Triangulation<2, 2>&, unsigned int, double, dealii::Triangulation<3, 3>&)’
   extrude_triangulation(const Triangulation<2, 2> &input,
   ^
/home/muhammad/installed_software/dealii-9.1.1/examples/Goal-oriented_mesh_adaptivity_in_elastoplasticity_problems/elastoplastic.cc:3959:30: error: invalid initialization of reference of type ‘dealii::Triangulation<3, 3>&’ from expression of type ‘dealii::parallel::distributed::Triangulation<2, 2>’
         extrude_triangulation(triangulation_2d,
                              ^
/home/muhammad/installed_software/dealii-9.1.1/examples/Goal-oriented_mesh_adaptivity_in_elastoplasticity_problems/elastoplastic.cc:107:3: note: in passing argument 4 of ‘void ElastoPlastic::extrude_triangulation(const dealii::Triangulation<2, 2>&, unsigned int, double, dealii::Triangulation<3, 3>&)’
   extrude_triangulation(const Triangulation<2, 2> &input,

Currently I am not changing anything in ".prm" and ".cc" and running it as it is downloading from the code gallery. May be I am missing some instructions to run it so any help or hint in this regard would be much appreciated. Thank you very much in advance!

[Wolfgang Bangerth]

Muhammad,
let me bring the original author of the program (Seyed Shahram Ghorashi) into
the loop as well.

> I ran the case of "Cantiliver_beam_3d" and it returned the results after
> couple of minutes as mentioned in the attachement (where per time step,
> seemingly only one iteration was performed in Newton method and in all
> increamental load steps the result values are zero for both stress and
> displacement).

Good question. I can confirm that that is what the program generates (though
at least the dual solution is non-zero), but as with all of the code gallery
programs, we just make them available for others to look at -- we don't
generally know what specifically they are doing, or whether they are correct.
You will have to ask Seyed Shahram Ghorashi this question.

> Secondly I tried to run the case "Thick_tube_internal_pressure" which first
> popped up the dialog message of :

> [...]

> Where as per my understanding of template based nature of deall.ii and also

> this case is 2d so if I change the *dim = 2* from originally *dim = 3* in main
> function of elastoplastic.cc

Yes, this is the correct solution.

> then the following error comes during compilation :

> [...]

Yes, I can confirm this as well. I will admit that I don't think this ever worked.

In the end, this is really the limitation of the code gallery. We see it as a
way for our users to share codes that may be of use to others in the
community, but we do not endorse them or verify the quality of these codes.
It's also quite possible that they *used* to work correctly with older deal.II
versions, but don't any more -- I could imagine that being the situation with
the first problem above.

Let's hope that Seyed has something to add that help you!

Best
WB

[Muhammad Mashhood]

Thank you Prof Wolfgang for your possible support and concern. I will be waiting for the further support from the author to successfully run it and to understand the workflow correctly.

[Muhammad Mashhood]

Dear Prof Wolfgang,

                              Hi! I am successful in running the code and it is without a doubt a nice addition to deal.ii code gallery. It reproduces the expected result for 3d Beam case now.

Since I did not find any thermomechanical coupled example so far in deal.ii therefore I want to contribute in the development and widening the scope of this code by coupling it with heat equation of step-26.

For this purpose , taking your lectures as a reference, Initially for testing purpose (i.e. before importing temperature solution from step-26) , I have added the thermal load in the cell_rhs as following:

              SymmetricTensor<2, dim> thermal_strain_qpoint;

              /////////////////////////// assigning thermal strain as +ve temperature change in whole body (i.e. all cells expanding equally) ////////////////////////

                                          for (unsigned int i=0; i<dim; ++i)
                                              for (unsigned int j=0; j<dim; ++j)
                                                  if (i==j)
                                                      thermal_strain_qpoint[i][j] = (1-0.3) * 600e-6;//1e-3 ;      // (1-nu) * alpha * del T
                                                  else
                                                      thermal_strain_qpoint[i][j] = 0*1e-3;  // no shear thermal strain in case of homogenous material properties

               ////////////////////////////////////////////////////////////////////////////////////////////////////////////

                  for (unsigned int j = 0; j < dofs_per_cell; ++j)
                    cell_matrix(i, j) += (stress_phi_i
                                          * fe_values[displacement].symmetric_gradient(j, q_point)
                                          * fe_values.JxW(q_point));

                  cell_rhs(i) += (
                                   ( stress_phi_i
                                     * incremental_strain_tensor[q_point] )
                                   -
                                   ( ( stress_strain_tensor
                                       * fe_values[displacement].symmetric_gradient(i, q_point))
                                     * tmp_strain_tensor_qpoint )
                                   +
                                   ( fe_values[displacement].value(i, q_point)
                                     * rhs_values_body_force )
                                   +
                                   ( fe_values[displacement].symmetric_gradient(i, q_point)
                                    * stress_strain_tensor * thermal_strain_qpoint   ) // adding the thermal loading force  produced from thermal strains

                                 ) * fe_values.JxW(q_point);

The solid geometry is in the form of rectangular plate (as shown in the attachment where mesh is even finer in the example being tested) with zero displacement (Dirchlet BC) at left edge and so far no mechanical traction load is applied.

In the result the displacement result is qualitatively as per the expectation i.e. the plate is expanding from other three ends but the System residuals in Newton iterations are of the order of 1e+4 or even more.

I do not have much experience in under concentration mathematical approach or I might be mistaking in writing the code for thermal loading force term as per deal.ii standards (because the rhs_values_body_force term already present in code is even bigger in numerical quantity than the thermal loading term which is added by me but still it lets system converge for small residuals of 1e-5 during solution) therefore I would be grateful to you for any guiding response in my this course of development. Thank you! 

[Muhammad Mashhood]

Similarly the example output log file is attached where the case is being simulated in which the gravity is acting for first two loading cases (which converge successfully) and from 3rd step the thermal load is started acting upon the plate and from this very loading step the Newton solution with residual: 115.471 is being accepted (and does not change for future Newton iterations) although the residuals for CG solver iterations and the difference of two consecutive solutions ( i.e. incremental displacements) is going even below 1e-6. 

[Wolfgang Bangerth]

Muhammad,

>                               Hi! I am successful in running the code and it
> is without a doubt a nice addition to deal.ii code gallery. It reproduces the
> expected result for 3d Beam case now.

Great. Can you say what you had to change to make things run? We should apply
these changes to the version that is available on the web so that others don't
have to figure this out like you had.

> Since I did not find any thermomechanical coupled example so far in deal.ii
> therefore I want to contribute in the development and widening the scope of
> this code by coupling it with heat equation of step-26.
> For this purpose , taking your lectures as a reference, Initially for testing
> purpose (i.e. before importing temperature solution from step-26) , I have

> added the thermal load in the /*cell_rhs*/ as following:
> [...]

> In the result the displacement result is qualitatively as per the expectation
> i.e. the plate is expanding from other three ends but the System residuals in
> Newton iterations are of the order of 1e+4 or even more.
> I do not have much experience in under concentration mathematical approach or
> I might be mistaking in writing the code for thermal loading force term as per

> deal.ii standards (because the /*rhs_values_body_force */term already present

> in code is even bigger in numerical quantity than the thermal loading term
> which is added by me but still it lets system converge for small residuals of
> 1e-5 during solution) therefore I would be grateful to you for any guiding
> response in my this course of development. Thank you!

I'm not quite sure what your question is, but for sure step-26 is not using
any physical units, so the solution there is not usable to define what the
right order of magnitude is.

In general, 1e+4 is not in itself a large number. You need to say what your
units are, and compare that with the typical size of what you are considering.
For example, if you consider the deformation of a bridge as you drive a truck
over it, and you choose to express the deformation in micrometers, then 1e4 is
not a large number. On the other hand, if you express the deformation in
lightyears, then 1e4 *is* a large number. So the question you need to ask
yourself is whether the residual is large *compared to what you would expect*.

Best
W.

[Muhammad Mashhood]

Dear Prof. Wolfgang,

                                Thanks for the guidance. I would be grateful to share the corrections we made for the others.

1) I concentrated upon running the 3d case of this code therefore the relevant changes made are in the following attached file. While main changes made by me and my colleague were to use set_boundary_id() rather than set_manifold_id  and some small modifications in refine_grid() function regarding distruibed_history_stress_field_comp and distruibed_history_strain_field_comp etc.

2) From the non linear residuals side after adding the thermal load in the rhs of elastoplastic problem, I am interested in values = 0 because here the equilibrium is being targeted in the body. So far I am also now able to get equilibrium from zero non linear residual values.

3) During my next step of combining the both codes (step-26) and (elastoplastic), I am facing some challenge that step-26 is so far in serial while the elastoplastic code is using the TrilinosWrapper::MPI::Vector. In your expert opinion what would be suitable to make them both communicate with each other to share the Vector<double> solution vector from step-26 to elastoplastic problem by running both at same imported mesh?

So far if I am correct, my idea is to convert step-26 problem also for parallel processing and then integrating it in the elastoplastic code's class i.e. both use the same triangulation. Any guiding methodology in this regard by sharing your experience or any suitable deal.ii example case etc. would be helpful.

Thank you in advance for your cooperation!