Crack propagation

[Seyed Ali Mohseni]

Dear Timo Heister and Thomas Wick,

I am trying to run your phase-field crack propagation example, namely the single edge notched tension test, with less than 9 Gauss points. 

Hence, the fe variable is initialized using FE_Q<dim>(degree) where degree cannot be 0. This means degree should always be >1 here in order to make things work, but how can one compute the examples with 1 or 4 Gauss points? Is the code not suited for what I am trying to achieve or where are the parts I have to change in order to make it work for less integration points. To be honest, I don't even see the point of having 9 integration points for a 4-noded quadrilateral or am I overseeing something?

Thank you in advance.

Kind regards,

S. A. Mohseni  

[Thomas Wick]

Dear S. A. Mohseni,

the Gauss points are assigned in that program as always in deal.II:

 QGauss<dim> quadrature_formula(degree + 2);

and in the argument, you can change to whatever number you want:

QGauss<dim> quadrature_formula(n);

with n=1 or n=2 or n = 3 for instance.

Best regards,

Thomas W.

-- 
++--------------------------------------------++
Dr. Thomas Wick
Maitre de conferences / Assistant Professor

Centre de Mathematiques Appliquees (CMAP)
Ecole Polytechnique
91128 Palaiseau cedex, France

Email:  thoma...@cmap.polytechnique.fr
Phone:  0033 1 69 33 4579
www:    http://www.cmap.polytechnique.fr/~wick/
++--------------------------------------------++
-- 

--
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.
For more options, visit https://groups.google.com/d/optout.

[Seyed Ali Mohseni]

Dear Thomas Wick,

I see. So I should change it manually everywhere. I thought the degree variable will do this automatically, but why did you then implement it like this? Is it for another purpose?

[Timo Heister]

> To be honest, I don't even see the point of having 9 integration points for
> a 4-noded quadrilateral or am I overseeing something?

Yes, you are missing some background information:
Quadrature formulas are exact for polynomials up to a certain degree.
For Gauss quadrature you need degree n to integrate polynomials of
degree 2n-1 exactly:
https://urldefense.proofpoint.com/v2/url?u=https-3A__en.wikipedia.org_wiki_Gaussian-5Fquadrature-23Fundamental-5Ftheorem&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=PDOwphRntO6TO2o4Lgd_VolY1Tdk0-uWL7LHVKKC9yY&s=GP1u_R42mj2e2VpImgPCnEHpalb6rChFLMnloN6wb_Q&e=

The integrals in the weak form involve products of shape functions of
a certain degree, so the quadrature should be of high enough degree to
integrate them with a small or no error. Maybe you can take a look at
https://urldefense.proofpoint.com/v2/url?u=http-3A__www.math.colostate.edu_-7Ebangerth_videos.676.33.5.html&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=PDOwphRntO6TO2o4Lgd_VolY1Tdk0-uWL7LHVKKC9yY&s=FGaRjjjHuJ-ybWRwN7eJRYI_p9YKmCjsIaS5mR2d458&e= or your
favorite FEM book.


--
Timo Heister
http://www.math.clemson.edu/~heister/

[Thomas Wick]

On 03/30/2017 03:38 PM, 'Seyed Ali Mohseni' via deal.II User Group wrote:

Dear Thomas Wick,

I see. So I should change it manually everywhere. I thought the degree variable will do this automatically, but why did you then implement it like this? Is it for another purpose?

It is just copy&paste from all the other deal.II-tutorials.

Just insert another global variable

unsigned int number_gauss_pt = XXX;

And then replace it once globally and for the future you again need only
to change one number in number_gauss_pt.

Best Thomas W.

Kind regards,

S. A. Mohseni

[Seyed Ali Mohseni]

Dear Thomas,

I tried what you suggested and it works. Thank you :)

There is a slight problem though when I try to compute the same example with the same boundary conditions for the unit_square_4.inp file. I am trying to run a one element test using your phase-field approach. It works, but there is no residual value and next to it I obtain a -nan below Reduction. Is it possible to run a one element test with the same settings or do I have to modify something (such as rhs computations or boundary value functions)?

Kind regards,

Seyed Ali

[Thomas Wick]

Dear Seyed Ali,

what do you mean by "one element test"?

Really to compute on one single element ?

This can never work:

1. First of all the phase-field will not be resolved
    because of the regularization parameter eps.

2. This contradicts any idea of numerical methods
    to have a reasonable number of mesh elements
    and then do mesh refinement with as many mesh elements as possible
    (w.r.t. to memory and computational cost)
    in order to approach the underlying continuous problem.

  
Therefore, the short answer is: no - a one element test is not possible
and it is also not reasonable to envisage it.

Best regards,

Thomas

-- 
++--------------------------------------------++
Dr. Thomas Wick
Maitre de conferences / Assistant Professor

Centre de Mathematiques Appliquees (CMAP)
Ecole Polytechnique
91128 Palaiseau cedex, France

Email:  thoma...@cmap.polytechnique.fr
Phone:  0033 1 69 33 4579
www:    http://www.cmap.polytechnique.fr/~wick/
++--------------------------------------------++
-- 

[Seyed Ali Mohseni]

Dear Thomas,

I would like to check the elastic energy you compute and compare it to my results. For a single element the computation in an elastic regime should be independent from the phase-field approach in primary steps.

BR,

Seyed Ali 

[Thomas Wick]

Dear Seyed Ali,

but this does not make sense since the finite
element approximation on one element is far from being meaningful.
Anything you obtain has a huge discretization error.
In addition, this result will not be independent of phase-field.

I would rather do it as we did (and usually done in numerics),
to compare the method to other published results or to
compute a solutions on a sequence of meshes.

For instance in the crack code, we compared our results
to Miehe et al. and Borden et al.

Best regards,

Thomas

-- 
++--------------------------------------------++
Dr. Thomas Wick
Maitre de conferences / Assistant Professor

Centre de Mathematiques Appliquees (CMAP)
Ecole Polytechnique
91128 Palaiseau cedex, France

Email:  thoma...@cmap.polytechnique.fr
Phone:  0033 1 69 33 4579
www:    http://www.cmap.polytechnique.fr/~wick/
++--------------------------------------------++
-- 

[Seyed Ali Mohseni]

Dear Thomas,

Exactly, I did that already. It seems it is identical, but I would like to see, if your energy variable psi_e is identical to mine. This is quite cumbersome for 500 elements or more to check. I mean the computation in the first steps is anyway independent from the phase-field value due to being in the elastic regime and the body has not been fractured yet. So, I assume it should not be a problem to just check the energy within the first steps. Of course you are right, the moment we reach the fracture toughness G_c cracks will start to grow and we have to separate between crack energy and bulk energy.

I cite your paper ;)

The results illustrate the expected behavior: as long as the crack does not grow (up to t = 0.0095 s), there is only an increase in bulk energy. Once the crack starts growing, bulk energy is dissipated into crack energy.

By the way Thomas, if you don't remember me, we already met after your GACM 2015 conference presentation and within the GAMM phase-field workshop :) Maybe you remember: Back then you already mentioned that h < epsilon to resolve the transition of the phase-field variable. Nevertheless, this is true for the cracked region, but if you have only elastic case the mesh size will not matter. This is just to allow a correct and accurate resolution and regularization of the phase-field crack region. 

[Thomas Wick]

Dear Seyed Ali,

now I understand better. What you could do is to disable
the phase-field functionality in the program or to
choose G_c very high such that the phase-field will never
be activated. Then you only have the elasticity part
and can compare.

By the way Thomas, if you don't remember me, we already met after your GACM 2015 conference presentation and within the GAMM phase-field workshop :)

Indeed. I did not make the relation ...

Maybe you remember: Back then you already mentioned that h < epsilon to resolve the transition of the phase-field variable. Nevertheless, this is true for the cracked region, but if you have only elastic case the mesh size will not matter. This is just to allow a correct and accurate resolution and regularization of the phase-field crack region.

Yes. This is correct. However, in case you do not know when the
crack starts propagating and phase-field values are going to zero
in the crack region, one has to be cautious with the number of
mesh cells (sufficiently enough), otherwise the relation h < epsilon
will be violated.

But again: even if you only have elasticity and thus do not
need to take care of h < epsilon; too coarse
mesh cells have a huge discretization error and are therefore
NOT accurate and thus questionable in their meaning!



Best regards,

Thomas

Kind regards,

S. A. Mohseni

[Seyed Ali Mohseni]

Dear Thomas,

Thank you. I checked the values and they are similar.

Within your predictor-corrector approach, I assume that you somehow exchange data between each step in order to be able to interpolate the old solution. Am I able to use the same functions to store energy data in each Gauss point within each element? Or am I forced to implement and use the quadrature point history approach from step-18? Due to the storage of all the values in each element, the computation becomes inefficient. It would be nice to be able to access data anywhere and anytime. Actually, if you are using Miehe's approach I am confident that you will have to store the energy somewhere in order to find the maximum value of the energy. Is that true?

Kind regards,

Seyed Ali 

[Thomas Wick]

Dear Seyed Ali,

On 04/01/2017 01:19 PM, 'Seyed Ali Mohseni' via deal.II User Group wrote:

Dear Thomas,

Thank you. I checked the values and they are similar.

Good.

Within your predictor-corrector approach, I assume that you somehow exchange data between each step in order to be able to interpolate the old solution.

Yes. But here we use the functionality that is already very well-developed in deal.II.

Am I able to use the same functions to store energy data in each Gauss point within each element? Or am I forced to implement and use the quadrature point history approach from step-18? Due to the storage of all the values in each element, the computation becomes inefficient. It would be nice to be able to access data anywhere and anytime. Actually, if you are using Miehe's approach I am confident that you will have to store the energy somewhere in order to find the maximum value of the energy. Is that true?

In principle you are right. If you use Miehe's approach with the strain-energy function to impose crack irreversibility, I would personally
use the idea from step-18. We however, inspired by Miehe's other paper developed our own approach
to impose the crack irreversibility constraint, namely by penalization. Either you really use penalization (simple,
augmented Lagrangian, etc. see my other papers with Wheeler/Mikelic and Lee etc.) or
you impose penalization by active set as down with Timo. But in all these approaches
we avoid indeed to store all energy values and create a history.


Best regards,

Thomas

Kind regards,

Seyed Ali 

[Seyed Ali Mohseni]

Dear Thomas,

I was wondering, why my computation is rather slow compared to your results within your paper. For instance regarding the shear test you mention that for global refinement the total time amounts to 5036 s. In my case each step takes around 200 s which would result in a total time of 30200 s (for 151 steps). What could be the reason?

Attached you find the terminal output from your approach.

[Thomas Wick]

Dear Seyed Ali,

no idea, but possible reasons :

1. Do you compute on the same mesh?

2. do you compile in release mode or debug mode?
    Release is faster.

3. Possibly you can optimize deal.Il with respect to the CPU etc.
    On my desktop computer at work, deal.II is 2 - 4 times slower than
    on my laptop for instance.

Best Thomas

[Seyed Ali Mohseni]

Dear Thomas,

1. Do you compute on the same mesh?

Yes, it is exactly the same mesh which is illustrated within your paper: 16384 cells (50115 DoFs). You can observe the same value also in the terminal output. I selected 4 global pre-refinement steps and 2 adaptive refinement cycles.

 

2. do you compile in release mode or debug mode?
    Release is faster.

Release mode.

 3. Possibly you can optimize deal.Il with respect to the CPU etc.

    On my desktop computer at work, deal.II is 2 - 4 times slower than
    on my laptop for instance.

I am using an Intel® Core™ i7-4770 CPU @ 3.40GHz × 8 processor. Hence, I assumed it would be even faster than your Intel(R) Core(TM) i5-3320M C...@2.60 GHz machine. What MPI solver did you use? Intel with OpenMPI? Since I am using MPICH maybe there are some differences. I mean the difference in speed is not magic and should be based on the processor type or the settings of the solver. Maybe I should check how fast my MPI solver works outside of deal.II in order to see where the problem lies. 

Best regards,

Seyed Ali  

[ha.ba...@gmail.com]

Hi

Each Core i7 CPU has maximum 4 physical cores. So, in MPI mode, solving the simulation using 4 core is faster than 8 cores in your computer.

 H.B

[Seyed Ali Mohseni]

Hi,

 

Each Core i7 CPU has maximum 4 physical cores. So, in MPI mode, solving the simulation using 4 core is faster than 8 cores in your computer.

This is not true. I checked it, it is slower.

Best,

Seyed Ali 

[Seyed Ali Mohseni]

Dear Thomas,

I found the problem. Your version works exactly like discussed in your paper, even faster with my Intel i7, namely around 1400 s. The timeconsuming part is the computation of my B-operator.

The B-operator computation is implemented inside assemble_system (bool residual_only):

for (; cell != endc; ++cell)

    if ( cell->is_locally_owned() )

    {

fe_values.reinit(cell);

local_matrix = 0;

local_rhs = 0;

int cell_index = cell->active_cell_index();

        ...

        // COMPUTE: B-Operator

        for (unsigned int q = 0; q < n_q_points; ++q) // loop over integration/quadrature points (GPs)

        {

             b_operators[cell_index][q] = Tensors::get_b_operator(fe_values, dofs_per_cell, q);

        }

        B = b_operators;

       
        // Old Newton iteration values

fe_values.get_function_values(rel_solution, old_solution_values);

fe_values.get_function_gradients(rel_solution, old_solution_grads);

        ...
     
    }

        

The variable B is a history variable I created to store the B-operator history and has been initialized in the main class function according to:

std::vector<std::vector<FullMatrix<double>>>B; // temporary cell storage solution (inefficient!)

The function get_b_operator() which is called within the Gauss points loop is implemented in Tensors:

template<int dim>

FullMatrix<double> get_b_operator (const FEValues<dim> &fe_values, const unsigned int dofs_per_cell, const unsigned int q)

{

FullMatrix<double> tmp(dim, GeometryInfo<dim>::vertices_per_cell);

// Remark: For vector-valued problems each column is a value for each value of the solution variable, hence here 3 DoFs per node (dofs_per_cell = 12)!

for (unsigned int i = 0; i < dofs_per_cell; i += 3)

{

const unsigned int index = i / 3;

// COMPUTE: B-operator (Remark: This version has to be extended for 3D!)

tmp[0][index] = fe_values.shape_grad_component(i, q, 0)[0];

tmp[1][index] = fe_values.shape_grad_component(i, q, 0)[1];

}

return tmp;

}

I assume (as I already wrote within my comments inside the above code) that this kind of cell storage for matrices turns out to be quite slow and inefficient. Jean-Paul offered me once to use the new CellDataStorage function, if I remember correctly. Isn't there a nice and elegant way to store data in each Gauss point and cell without loss of computational time? Or am I able to pass the B-operator to another function somehow after I have solved the system? Something similar to your predictor-corrector scheme. It seems you somehow pass data more efficiently from one function to another function even after you solved the system. Do you also pass cell data from a function before solving to a function after solving?

Best,

Seyed Ali

[Thomas Wick]

Dear Seyed Ali,

very good!

However, when you write that you solve "exactly" my/our (Timo and I)
code then I would assume it is really "exactly" that version.

When you add other things, of course this might accelerate or slow
down the code :)

Next time, please write first if you solve "exactly" the same code or not :)

Anyhow: good that you found the reason.

Best regards,

Thomas

-- 
++--------------------------------------------++
Dr. Thomas Wick
Maitre de conferences / Assistant Professor

Centre de Mathematiques Appliquees (CMAP)
Ecole Polytechnique
91128 Palaiseau cedex, France

Email:  thoma...@cmap.polytechnique.fr
Phone:  0033 1 69 33 4579
www:    http://www.cmap.polytechnique.fr/~wick/
++--------------------------------------------++
-- 

[Seyed Ali Mohseni]

Dear Thomas,

Sure, my bad. I am sorry for the inconvenience it may have caused you and Timo. Unfortunately, I was just paying attention to the input text file and oversaw the changes I actually made. Won't happen again. ;)

Regarding the speedup, I just figured out how to make use of the CellDataStorage function which has been recently implemented within deal.II. First of all, I want to thank Denis Davydov and Jean-Paul Pelteret for implementing it so nice and stable. :) Without their recent work, I probably would have to find a way to write it myself somehow. Regarding the CellDataStorage approach: Is it possible to store a FullMatrix<double> within the MyQData struct function? If yes, do I have to write separate functions for pack_values and unpack_values functions?

Kind regards,

Seyed Ali

[Denis Davydov]

Hi Seyed,

On Sunday, April 9, 2017 at 8:18:10 PM UTC+2, Seyed Ali Mohseni wrote:

Dear Thomas,

Sure, my bad. I am sorry for the inconvenience it may have caused you and Timo. Unfortunately, I was just paying attention to the input text file and oversaw the changes I actually made. Won't happen again. ;)

Regarding the speedup, I just figured out how to make use of the CellDataStorage function which has been recently implemented within deal.II. First of all, I want to thank Denis Davydov and Jean-Paul Pelteret for implementing it so nice and stable. :) Without their recent work, I probably would have to find a way to write it myself somehow. Regarding

thanks for your kind words. I am glad that our class is useful to others. 

 

the CellDataStorage approach: Is it possible to store a FullMatrix<double> within the MyQData struct function? If yes, do I have to write separate functions for pack_values and unpack_values functions?

you only need those if you want to transfer quadrature data during refinement. But yes, if you store a FullMatrix<> at quadrature data, then you would need to 

implement those functions to enable transfer of the data during mesh refinement. All in all it's not different from storing several scalar valued data at each quadrature point.

Projection/interpolation of data is done independently for each scalar anyway, regardless of its interpretation (scalar, component of a vector, component of a second order tensor, etc). 

Note that this may not be what you want in the context of crack propagation because what we do to transfer is assume that the quadrature data is continuous across the element (thus the name ContinuousQuadratureDataTransfer).

Regards,

Denis.

[Seyed Ali Mohseni]

Dear Denis Davydov,

I am trying to implement the pack_values() and unpack_values() functions for a FullMatrix<double>, but the function itself in quadrature_point_data.h is defined using a vector of double values only. Nevertheless, my initial approach would be the following:

struct MyQData: public TransferableQuadraturePointData

{

double psi;

FullMatrix<double> B;

unsigned int number_of_values () const

{

return 2;

}

virtual void pack_values (std::vector<double> &scalars) const

{

Assert(scalars.size()==1, ExcInternalError());

scalars[0] = psi;

}

virtual void pack_values (std::vector<FullMatrix<double>> &matrices) const

{

Assert(scalars.size()==1, ExcInternalError());

matrices[0] = B;

}

virtual void unpack_values (const std::vector<double> &scalars)

{

Assert(scalars.size() ==1, ExcInternalError());

psi = scalars[0];

}

virtual void unpack_values (const std::vector<FullMatrix<double>> &matrices)

{

Assert(matrices.size() ==1, ExcInternalError());

B = matrices[0];

}

};

Is the above correct? Will it work during refinement? I checked it without refinement and it seems to give me correct values at each Gauss point after data storage.

Note that this may not be what you want in the context of crack propagation because what we do to transfer is assume that the quadrature data is continuous across the element (thus the name ContinuousQuadratureDataTransfer).

Yes, but what choice do I have? Is there another computationally efficient approach? 

This means probably that the PointHistory function data storage based on step-18 has the same drawback for discontinuous problems, am I right?  

I would be gladly willing to implement the DiscontinuousQuadratureDataTransfer function, if it is somehow possible to accomplish in deal.II. I am pretty sure it will be beneficial for others also in many ways. ;) 

[Denis Davydov]

That’s not the right way to do it, you should have something along the lines:

unsigned int number_of_values () const

{  

    return 1 + <whatever_is_the_number_of_elements_in_a_matrix>;

}

virtual void pack_values(std::vector<double> &values) const

{

    values[0] = psi;

    values[1] = B[0][0];

    values[2] = B[0][1];

    ...

}

and same for unpack but other way around (setting internal data from an input vector).

The point is that we can not guess what type of data users are storing. But we don’t really need to.

We only need to know how many independent scalar fields would the data correspond do.

That’s what is returned by number_of_values(). The same resasons pack_values() are NOT called pack_scalars() etc

Regards,

Denis.

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/rgv4YGx5VJE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.

[Jack]

Hi S. A. Mohseni,

I'm also interested on the cracking mechanics, where is such an example for crack propagation ? 

Thanks !

在 2017年3月30日星期四 UTC+8下午9:13:16，Seyed Ali Mohseni写道：

[Daniel Arndt]

Jack,

have a look at https://github.com/tjhei/cracks

Best,
Daniel

[Jack]

Hi Daniel,

Thank you very much!

Regards,

Jack

在 2017年5月14日星期日 UTC+8下午5:07:26，Daniel Arndt写道：

[fea...@gmail.com]

Timo,

If I use your user element, I have to use OPEN MPI? now I have some issues with Open MPI in Deal.ii

Regards,
Yaakov

[Timo Heister]

> If I use your user element, I have to use OPEN MPI? now I have some issues
> with Open MPI in Deal.ii

What do you mean by "user element"? The example code in question
requires deal.II to be configured with MPI. What vendor you use
(OpenMPI, MPICH, ...) is up to you.

--
Timo Heister
http://www.math.clemson.edu/~heister/

[fea...@gmail.com]

Dear Prof. Heister,

I would like to just test  isotrope formulation of phase -field model (no compression/tension modification), how could I modify the codes (in a simple way)?

I am sorry that I am just a beginner of DealII.

Best

[Thomas Wick]

Hello,

there are two flags in the parameter file that you need to change:

set Decompose stress in rhs               = 0.0
set Decompose stress in matrix           = 0.0


Thomas W.

-- 
++--------------------------------------------++
Prof. Dr. Thomas Wick
Institut für Angewandte Mathematik (IfAM)
Leibniz Universität Hannover
Welfengarten 1
30167 Hannover, Germany

Tel.:   +49 511 762 3360
Email:  thoma...@ifam.uni-hannover.de
www:    http://www.ifam.uni-hannover.de/wick

www:    http://www.cmap.polytechnique.fr/~wick/
++--------------------------------------------++
-- 

--

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.

[fea...@gmail.com]

Dear Prof. Wick,

thanks a lot for your support.

I will learn how these commands work in the proposed phase-field user element.

Kind regards,
Yaakov

[fea...@gmail.com]

Dear Prof. Wick,

I have used isotrope formulation for miehe shear loading (without local refinement). I ca not see two crack branching which is described in the article (A review on phase-field models of brittle fractureand a new fast hybrid formulation)

I attach the test results.

Thanks for your answer!

Kind regards,
Yaakov

On Saturday, February 3, 2018 at 7:20:27 AM UTC+1, Thomas Wick wrote:

[Thomas Wick]

Dear Yaakov,

which article do you mean? Please give the exact reference including author names.

I do not know a priori whether they have different material parameters, another
stress-strain splitting, etc.

The reason for different results can be anything. One needs to do a careful 1-by-1
comparison.

Best regards,

Thomas W.

[Thomas Wick]

Dear Yaakov,

which article do you mean? Please give the exact reference including author names.

I do not know a priori whether they have different material parameters, another
stress-strain splitting, etc.

The reason for different results can be anything. One needs to do a careful 1-by-1
comparison.

Best regards,

Thomas W.