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