Problem with postprocessor

[Seyed Ali Mohseni]

Dear Thomas Wick, Dear Timo Heister,

I wrote an additional postprocessor in your existing phase-field code to allow postprocessing of strain, stress or elastic energy. Unfortunately, it seems like the STRAIN_XX and STRAIN_XY is not increasing in each step while the STRAIN_YY increases correctly.

The numerical example is a simple 2d cube of length and width 1.0 consisting of 1 linear element only. The material properties are the same from the Miehe tension example. The specimen is loaded exactly the same way as in the tension experiment.

Of course, I am aware of the fact that the element size is totally wrong for a phase-field simulation. I am merely trying to check the strain in a purely elastic case before crack initiation.   

Is it due to the predictor-corrector nature or something which causes this problem?

In my solid_mechanics code written in deal.II it works correctly, the following postprocessor implementation:

template<int dim>

class Postprocessor: public DataPostprocessor<dim>

{

public:

Postprocessor ();

void compute_derived_quantities_vector (const std::vector<Vector<double> > &uh, const std::vector<std::vector<Tensor<1, dim> > > &duh, const std::vector<std::vector<Tensor<2, dim> > > &dduh, const std::vector<

Point<dim> > &normals, const std::vector<Point<dim> > &evaluation_points, std::vector<Vector<double> > &computed_quantities) const;

virtual std::vector<std::string> get_names () const;

virtual std::vector<DataComponentInterpretation::DataComponentInterpretation> get_data_component_interpretation () const;

virtual UpdateFlags get_needed_update_flags () const;

private:

void print_tensor (const Tensor<2, dim>& tensor, char* name) const;

};

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

template<int dim>

Postprocessor<dim>::Postprocessor ()

{

}

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

template<int dim>

void Postprocessor<dim>::compute_derived_quantities_vector (const std::vector<Vector<double> > &uh, const std::vector<std::vector<Tensor<1, dim> > > &duh, const std::vector<

std::vector<Tensor<2, dim> > > &/*dduh*/, const std::vector<Point<dim> > &/*normals*/, const std::vector<Point<dim> > &evaluation_points, std::vector<

Vector<double> > &computed_quantities) const

{

//TODO: Postprocessing has not been optimized for 3D yet.

// Number of quadrature points (interior)

const unsigned int n_q_points = uh.size();

// CHECK: Evaluation points

const std::vector<Point<dim> > EP = evaluation_points;

// std::cout << "\nEVALUATION POINTS\n" << std::endl;

// for (unsigned int i = 0; i < n_q_points; ++i)

// {

// for (unsigned int j = 0; j < dim; ++j)

// std::cout << "  " << EP[i][j] << " ";

// std::cout << std::endl;

// }

// std::cout << std::endl;

// Constitutive matrix

SymmetricTensor<4, dim> C = Tensors::get_elasticity_tensor<dim>(FractureMechanics<dim>::public_lambda, FractureMechanics<dim>::public_mu);

for (unsigned int q = 0; q < n_q_points; ++q)

{

std::cout << "\n - in EVALUATION point " << q + 1 << " - \n" << std::endl;

Tensor<2, dim> grad_u;

SymmetricTensor<2, dim> eps, sigma;

double eps_ii, eps_ij, sigma_ii, sigma_ij;

// ===========================[ STRAINS ]============================

for (unsigned int i = 0; i < dim; ++i)

{

int j = i + 1;

grad_u[i] = duh[q][i];

std::cout << "DUH 0: " << duh[q][0] << std::endl;

std::cout << "DUH 1: " << duh[q][1] << std::endl;

eps = 0.5 * (grad_u + transpose(grad_u));

eps_ii = eps[i][i];

if ( j < dim )

eps_ij = eps[i][j];

// eps_ij = 2.0 * eps[i][j];

// std::cout << "STRAIN " << i << i << ": " << eps_ii << std::endl;

// std::cout << "STRAIN " << i << j << ": " << eps_ij << std::endl;

computed_quantities[q](i) = eps_ii;

computed_quantities[q](j) = eps_ij;

}

// ==========================[ STRESSES ]============================

for (unsigned int i = 0; i < dim; ++i)

{

int j = i + 1;

sigma = C * eps;

sigma_ii = sigma[i][i];

if ( j < dim )

sigma_ij = sigma[i][j];

computed_quantities[q](dim + i + 1) = sigma_ii;

computed_quantities[q](dim + j + 1) = sigma_ij;

}

// =======================[ ELASTIC ENERGY ]=========================

// double psi = 0.5 * FractureMechanics<dim>::public_lambda * trace(eps) * trace(eps) + FractureMechanics<dim>::public_mu * eps * eps;

double psi = 0.5 * (eps[0][0] * sigma[0][0] + eps[1][1] * sigma[1][1] + 2.0 * eps[0][1] * sigma[0][1]);

// double psi = 0.5 * scalar_product(sigma,eps);

// std::cout << std::endl << "DISPLACEMENT GRADIENT" << std::endl;

// for (unsigned int i = 0; i < dim; ++i)

// {

// for (unsigned int j = 0; j < dim; ++j)

// std::cout << " " << std::setprecision(6) << std::fixed << grad_u[i][j] << "  ";

// std::cout << std::endl;

// }

// std::cout << std::endl;

// print_tensor(eps, "STRAIN TENSOR");

// print_tensor(sigma, "STRESS TENSOR");

computed_quantities[q](dim * 2 + 2) = psi;

}

}

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

template<int dim>

std::vector<std::string> Postprocessor<dim>::get_names () const

{

std::vector<std::string> names;

std::string xyz = "xyz";

// =================[ STRAIN NAMES ]==================

// Diagonal strain components

for (int i = 0; i < dim; i++)

{

std::stringstream ss;

ss << "strain_" << xyz[i] << xyz[i];

std::string s = ss.str();

names.push_back(s);

}

// Symmetric strain components (3d may not work!)

for (int i = 0, j = 1; j < dim; j++)

{

std::stringstream ss;

ss << "strain_" << xyz[i] << xyz[j];

std::string s = ss.str();

names.push_back(s);

}

// =================[ STRESS NAMES ]==================

// Diagonal stress components

for (int i = 0; i < dim; i++)

{

std::stringstream ss;

ss << "stress_" << xyz[i] << xyz[i];

std::string s = ss.str();

names.push_back(s);

}

// Symmetric stress components (3d may not work!)

for (int i = 0, j = 1; j < dim; j++)

{

std::stringstream ss;

ss << "stress_" << xyz[i] << xyz[j];

std::string s = ss.str();

names.push_back(s);

}

// ==================[ ENERGY NAME ]==================

names.push_back("elastic_energy");

return names;

}

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

template<int dim>

std::vector<DataComponentInterpretation::DataComponentInterpretation> Postprocessor<dim>::get_data_component_interpretation () const

{

std::vector<DataComponentInterpretation::DataComponentInterpretation> interpretation(dim, DataComponentInterpretation::component_is_scalar);

interpretation.push_back(DataComponentInterpretation::component_is_scalar);

interpretation.push_back(DataComponentInterpretation::component_is_scalar);

return interpretation;

}

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

template<int dim>

UpdateFlags Postprocessor<dim>::get_needed_update_flags () const

{

return update_values | update_gradients | update_quadrature_points;

}

Also I wonder why the elastic energy is increasing to a large value such as 4000 or more whereby it should be around 200. But that is another problem which only occurs in a state where the crack propagates. Hence, when I compare the energies or strains from another code, in the elastic regime they are almost similar, but when the crack starts to propagate they differ. 

Hope you have an idea what could be causing such observations.

Kind regards,

S. A. Mohseni

[Thomas Wick]

Dear S. A. Mohseni,

On 05/31/2017 01:11 PM, 'Seyed Ali Mohseni' via deal.II User Group wrote:

Dear Thomas Wick, Dear Timo Heister,

I wrote an additional postprocessor in your existing phase-field code to allow postprocessing of strain, stress or elastic energy. Unfortunately, it seems like the STRAIN_XX and STRAIN_XY is not increasing in each step while the STRAIN_YY increases correctly.

This is difficult to say and heavily depends on the loading conditions. The stress has different components, but not all
of them increase for a specific loading.

The numerical example is a simple 2d cube of length and width 1.0 consisting of 1 linear element only. The material properties are the same from the Miehe tension example. The specimen is loaded exactly the same way as in the tension experiment.

Of course, I am aware of the fact that the element size is totally wrong for a phase-field simulation. I am merely trying to check the strain in a purely elastic case before crack initiation.  

I stay with my previous opinion some weeks ago: it does not make sense to run a FE simulation only on one element !!!
Where is the problem to have 100 or 1000 elements? This simulation (see step-3 for Laplace) also would run
only for seconds.

Is it due to the predictor-corrector nature or something which causes this problem?

This is easy to check: just disable predictor-corrector and run with uniform mesh refinement and then you see
whether predictor-corrector causes the problem.

But anyhow: when you work only on one element - there is no predictor-corrector because otherwise you would have > 1 elements.



Best regards,

Thomas W.

--
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.

-- 
++--------------------------------------------++
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 already run the example with uniform mesh, hence global refinement with 0 refinement cycles.

The problem is, a specimen with 100 or 1000 elements is difficult to check due to the terminal output being flooded.

I could just output the result for one element, but then wheres the difference. 

Don't you ever do some simple FE patch examples?

[Thomas Wick]

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

Dear Thomas,

I already run the example with uniform mesh, hence global refinement with 0 refinement cycles.

The problem is, a specimen with 100 or 1000 elements is difficult to check due to the terminal output being flooded.

This I understand, but you could write everything into a file ...

I could just output the result for one element, but then wheres the difference.

See below.

Don't you ever do some simple FE patch examples?

Never. Because simply the theory of finite elements tell you that the discretization error is so huge
that any result is nearly meaningless.


If I am interested in stress values or the stress tensor, but I am not sure if this
will help you or if you are interested, is to compute typical quantities of interest.
That is not only a point-wise stress, but for example a line integral:

\int_{part of boundary} \sigma\cdot n \ds.

Or indeed you compute the stress in a specific point in the domain
but really using more than one element.

Many examples and benchmarks for elasticity and plasticity and also quantities
of interest (goal functionals) are given in the book:

http://eu.wiley.com/WileyCDA/WileyTitle/productCd-0471496502.html


For instance two examples from this book with point values of stresses (e.g., \sigma_yy),
again on more than one element, are provided in the deal-DOpElib library, see pages 42 - 45 on

http://wwwopt.mathematik.tu-darmstadt.de/dopelib/description_full.pdf


Best Thomas

Kind regards,

S. A. Mohseni

--
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.

[Tobi Young]

I'm going to jump in with one of my random comments uninvited. Hopefully you don't mind. :-)

I already run the example with uniform mesh, hence global refinement with 0 refinement cycles.

The problem is, a specimen with 100 or 1000 elements is difficult to check due to the terminal output being flooded.

I think what Thomas is kindly trying to point out is, that a single cell is a useless test case for almost all possible test cases that deal with practical problems. 

It seems your case is one of many examples.

Try 16 cells. Look at the data and then try 32 cells and compare. Maybe write a procedure that will look at the data for you and put useful results into a file for you to look at? Can you plot results with gnuplot (for example)?

Machines are stupid, they only can do what you tell them to do! So, why not reserve a day or so to calmly write the algorithms needed to extract the data you want in a way you can visualise?

You need to do your analysis in some way. You can not expect a machine to do things for you. You have to instruct her what it is you want to be done. Write the code...   :-)

That is scientific computing. 

I could just output the result for one element, but then wheres the difference. 

Big difference in numerical analysis of this kind! ;-)

I am not being unkind - though, maybe its time to get dirty, write some code, and look at the numbers, in order to figure out where the problem lies. If you can do that, there are alot of dealii users and developers that will help you out, and you'll get there in the end. :-)

Best,

   Toby

[Thomas Wick]

Hi Toby,

I appreciate your additional comments.

Thanks and best,

Thomas

--
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]

Hi again,

Nice words, indeed. Thank you. I will follow your advice :)

But to get back to the topic: Unfortunately, you misunderstood my problem. The issue is with the displacement gradient values stored in "duh" which is computed by deal.II and given as input for the postprocessing tasks to be written.

I compared the output of solution gradients "duh" for both codes (my own solid_mechanics code and Thomas' code) using the same implementation for the postprocessor. 

The fascinating thing is that my code works and it gives correct results in my benchmark since I validated it with another FE program of our own to be 100 % sure.

But the code written by Thomas and Timo somehow doesn't give correct duh values.

In the first steps both codes agree, but the duh values won't increase according to the loading, only the "yy" component in loading direction increases.

This is strange.

I am not completely familiar with Thomas' code that's why I posted here. Otherwise debugging my own code is what I do daily.

I hope you could explain me what could cause such a behavior in your implementation that deal.II gives different values for duh although the same postprocessor implementation is being used.

Is it because of the coupled formulation?

Thank you for your help so far.

Kind regards,

Seyed Ali 

[Thomas Wick]

I hope you could explain me what could cause such a behavior in your implementation that deal.II gives different values for duh although the same postprocessor implementation is being used.

Is it because of the coupled formulation?

I really would think so. What you could do is to disable the phase-field variable
such that you have in both codes really only elasticity and have a fair comparison.

How do you do this?

You could erase in our code (Heister/Wick) all the phase-field (pf) appearances in the solid
mechanics equations. Then the phase-field values do not enter any more
and led to wrong results with the fracture.


Best Thomas

Thank you for your help so far.

Kind regards,

Seyed Ali 

--
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]

I did another trick: Increasing G_c has the same effect to achieve purely elastic behavior. So, I chose a high enough G_c value and it works!

I obtain the same results and the duh values grow correctly. 

Still I cannot understand fully why when cracking is initiated, there is no increase of the strains except the strain in the loaded direction.

Maybe you are right, I try to compare a bigger example, but I already compared two phase-field codes.

I have my own version with staggered scheme and it was not identical to your results after crack initiation.

As mentioned before, in the elastic regime everything is identical and fine.

The only thing I changed a bit is the integration from quadratic to linear, but I am not sure, if seeting quadrature_formula(2) alone is enough for linear shape functions. You also had face_quadrature_formula or this Lobatto integration.

Would you be so kind and explain what lines I have to change in order to obtain linear shape functions and not quadratic like set in the original version?

How many changes and the code would help a lot.

Just to be sure I did it correct.

Thank you.

Best,

Seyed Ali 

[Thomas Wick]

On 05/31/2017 05:00 PM, 'Seyed Ali Mohseni' via deal.II User Group wrote:

I did another trick: Increasing G_c has the same effect to achieve purely elastic behavior. So, I chose a high enough G_c value and it works!

I obtain the same results and the duh values grow correctly.

Excellent. This is indeed also possible. Very good idea.

Still I cannot understand fully why when cracking is initiated, there is no increase of the strains except the strain in the loaded direction.

Maybe you are right, I try to compare a bigger example, but I already compared two phase-field codes.

I have my own version with staggered scheme and it was not identical to your results after crack initiation.

As mentioned before, in the elastic regime everything is identical and fine.

I made the experience that there can be large differences between staggered, quasi-monolithic (Heister/Wick) and
fully monolithic (I have some recent studies on this).

Moreover, sometimes it makes a difference how you impose the irreversibility condition: via strain history or
penalization.

The only thing I changed a bit is the integration from quadratic to linear, but I am not sure, if seeting quadrature_formula(2) alone is enough for linear shape functions. You also had face_quadrature_formula or this Lobatto integration.

Would you be so kind and explain what lines I have to change in order to obtain linear shape functions and not quadratic like set in the original version?

You need to change the number in FE_Q( ....) for another degree. But in our code,
it should be already linear if I remember well.


Best regards,

Thomas

How many changes and the code would help a lot.

Just to be sure I did it correct.

Thank you.

Best,

Seyed Ali 

--
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.

[Thomas Wick]

On 05/31/2017 05:00 PM, 'Seyed Ali Mohseni' via deal.II User Group wrote:

The only thing I changed a bit is the integration from quadratic to linear, but I am not sure, if seeting quadrature_formula(2)

In

https://github.com/tjhei/cracks/blob/master/cracks.cc

line 982,

you need to change the "degree", but it should be already "1".

And the quadrature formula - what you mentioned above - should be sufficiently high, "2" should do the job.


Best Thomas

alone is enough for linear shape functions. You also had face_quadrature_formula or this Lobatto integration.

Would you be so kind and explain what lines I have to change in order to obtain linear shape functions and not quadratic like set in the original version?

How many changes and the code would help a lot.

Just to be sure I did it correct.

Thank you.

Best,

Seyed Ali 

--
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,

Would you be so kind to tell me how am I able to create the load-displacement curves like in your paper? Did you use gnuplot?

And does deal.II have a built-in function that allows one to output the necessary data directly to gnuplot or to a textfile?

Thank you in advance.

Kind regards,

Seyed Ali

[Thomas Wick]

Hi Seyed,

I am old-fashioned in this sense: indeed yes, I just plot the output into a file:

make run | tee my_results.txt

Then I `grep' the quantity (for example PStress) that I want and then plot with gnuplot.

Of course it would be a bit more efficient to write all output of interest directly into
a file. Honestly this is not so difficult either with deal.II functions or directly
C++ (fstream, write, etc.).

In step-7 you will get an idea how to use deal.II functions for latex output.

Best Thomas

--
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.