Averaging derived solution quantities (mises stress)

[Konrad Schneider]

Dear all,

I am trying to generate a mises stress contour plot out of an elasticity simulation.
So far I managed to postprocess derived quantities from the displacement solution like displacement gradients to get strains and eventually stress for computation of the mises stress. However, when I view the contour plot of the von mises stress, the picture looks rather ugly (see attached Mises_quadrature_points.png).

with the discontinuities between the cells clearly visible. If I understand it correct, dealii gives me the gradient values of the solution at quadrature points. When I compute the mises stress and pass it to a DataOut-object it attaches these values to the vertices of each cell. What happens if there are more quadrature points than vertices in a cell? How are they then mapped to the vertices.
To polish the contour plot I tried to mimic the procedure of Step-18, where the stress is averaged in each cell. To this end I used the following piece of code:

Vector<double> mises_avg(this->triangulation.n_active_cells());

{

QGauss<dim> quadrature_formula(this->fe.degree + 1);

FEValues<dim> fe_values(this->mapping,

this->fe,

quadrature_formula,

update_gradients);

const unsigned int n_q_points = quadrature_formula.size();

std::vector<std::vector<Tensor<1, dim>>> u_gradients( n_q_points, std::vector<Tensor<1, dim>>(dim));

const SymmetricTensor<4, 3> C_0=get_stress_strain_tensor<3>(this->E,this->nu);

for (auto &cell : dof_handler.active_cell_iterators()){

fe_values.reinit(cell);

fe_values.get_function_gradients(u_glob, u_gradients);

double avg_mises_stress = 0.;

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

Tensor<2, dim> grad_u;

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

grad_u[d] = u_gradients[q][d];

const SymmetricTensor<2, dim> strain_dim = symmetrize(grad_u);

SymmetricTensor<2, 3> strain;

if (dim<3){// plain strain case

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

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

strain[k][l] = strain_dim[k][l];

}

const SymmetricTensor<2, 3> stress= C_0 * strain;

avg_mises_stress += std::sqrt(3./2.)*deviator(stress).norm();

}

mises_avg(cell->active_cell_index()) = avg_mises_stress/n_q_points;

}

}

data_out.add_data_vector(mises_avg, "mises_avg");

To receive the following picture:

Not much of an improvement.

Is there a possibility to get a smoother output? I thought of averaging the individual cell contributions connected to a vertex.

I would be grateful for any advice.

[Wolfgang Bangerth]

Simon,

> So far I managed to postprocess derived quantities from the displacement
> solution like displacement gradients to get strains and eventually stress
> for computation of the mises stress. However, when I view the contour plot
> of the von mises stress, the picture looks rather ugly (see attached

> Mises_quadrature_points.png).Mises_quadrature_points.png with the

> discontinuities between the cells clearly visible. If I understand it
> correct, dealii gives me the gradient values of the solution at quadrature
> points. When I compute the mises stress and pass it to a DataOut-object
> it attaches these values to the vertices of each cell. What happens if
> there are more quadrature points than vertices in a cell? How are they then
> mapped to the vertices.

How do you create this output? Via a DataPostprocessor object? In that case,
the gradients are computed in the vertices of cells, not quadrature points.

To me, the plot does not actually look particularly badly. Your mesh is just
too coarse :)

> To polish the contour plot I tried to mimic the procedure of Step-18, where
> the stress is averaged in each cell. To this end I used the following

> piece of code: [...]
> To receive the following picture: Screenshot from 2023-03-15 00-45-46.png

> Not much of an improvement. Is there a possibility to get a smoother
> output? I thought of averaging the individual cell contributions connected
> to a vertex.

You can do that, and I believe that there is a function somewhere in either
namespace DoFTools or VectorTools that averages cell-based quantities to
node-based quantities.

In reality, though, if you *really* want a smooth field, then you want to use
the Zienkiewicz-Zhu "recovery" procedure. It is a particular kind of averaging
that produces a stress (gradient) field that can be shown to be
one-order-more-accurate than one would naively expect. Of course, the
alternative is to just do what you are already doing but on a finer mesh :-)

Best
W.

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

[Konrad Schneider]

Dear Wolfgang,

yes the mesh is quite coarse, but I want to investigate the functionality of how deal.ii outputs stresses and how these stresses are then displayed via for instance paraview. To really see the effects I intentionally use the coarse mesh.
So far, I actually do create a (derived) DataPostprocessor-Object (see the thinned out and altered step-8 tutorial attached to this message) similar to the procedure in step-32. You wrote, that the gradients are computed at the vertices. How so? If I look at line 291 of my code or at step-32, there is a loop over all quadrature points (  for (unsigned int q = 0; q < n_quadrature_points; ++q) ) which gets called for every cell. In the scope of this loop the shape function gradients at the quadrature points are used to compute the desired output variable by a special calculation rule. However, these are not the vertices. If the gradients are somehow again computed at the vertices, which might be an internal procedure of the DataPostprocessor-functionality of dealii, how would this procedure know what the special calculation rule of the desired output would look like?

Concerning your second note, thanks for pointing me in the direction of  the Zienkiewicz-Zhu recovery procedure. The procedure is straightforward to me, however its implementation in deal.ii not (yet). There are some functions like the VectorTools::project() function, which is supposed to give a L2-projections. 
However, I am not quite sure how to use them. Therefore I tried to implement the  Zienkiewicz-Zhu recovery procedure myself. But there is a fundamental problem I could not solve so far. In elasticity problems for 2D or 3D we have a vector-valued solution. Therefore, in 2D for a triangulation with n Nodes we have 2n DoFs. If I want to add the projected scalar function (e.g. the mises stress) in form of a command like

data_out.add_data_vector(mises_stress_recovery, "mises_recovery");

to the output, there is the problem that the mises_stress_recovery vector only has a dimension of n, where as data_out wants a vector of dimension 2n. How to deal with that issue?
I tried to generate a new dof_handler via 

DoFHandler<dim> dof_handler_recovery(this->triangulation);

FESystem<dim> fe_recovery(FE_Q<dim>(this->fe_degree),1); // here dim=1 since we only have on dof per node

MappingQ<dim> mapping_recovery(this->fe_degree);

AffineConstraints<double> constraints_recovery;

dof_handler_recovery.distribute_dofs(fe_recovery);

and then assemble a mass matrix and a rhs to solve for the mises_stress_recovery vector. However, when looping over the cells, I am not sure how to reference to the cells. Should I use the cells of the newly created DoFHandler via 

for (const auto &cell : dof_handler_recovery.active_cell_iterators())

or via the DoFHandler of the original problem solution via

for (const auto &cell : dof_handler.active_cell_iterators()).

Do I really have to create a new DoFHandler? If not, how would one assemble the mass matrix, since for that we need the an local_dof_indices object like std::vector<types::global_dof_index> that points only to every second dof index or so.
I would be grateful for some tips or directions.

Best 
Konrad

[Konrad Schneider]

there was a small bug in my code I posted above. the following code is running fine.

[Wolfgang Bangerth]

Konrad,

> yes the mesh is quite coarse, but I want to investigate the functionality of
> how deal.ii outputs stresses and how these stresses are then displayed via for
> instance paraview. To really see the effects I intentionally use the coarse mesh.
> So far, I actually do create a (derived) DataPostprocessor-Object (see the
> thinned out and altered step-8 tutorial attached to this message) similar to
> the procedure in step-32. You wrote, that the gradients are computed at the
> vertices. How so? If I look at line 291 of my code or at step-32, there is a
> loop over all quadrature points (

> for(unsignedintq=0;q<n_quadrature_points;++q) ) which gets called for every
> cell.

Well, there is a loop over all of the evaluation points for which DataOut
gives your function the solution. In that function, you *call* these
quadrature points, but they are just some evaluation points, unrelated to what
you use when you compute integrals.

> In the scope of this loop the shape function gradients at the quadrature
> points are used to compute the desired output variable by a special
> calculation rule. However, these are not the vertices. If the gradients are
> somehow again computed at the vertices, which might be an internal procedure
> of the DataPostprocessor-functionality of dealii, how would this procedure
> know what the special calculation rule of the desired output would look like?

Can you clarify this question? What are the "special rules"?

> Concerning your second note, thanks for pointing me in the direction of  the
> Zienkiewicz-Zhu recovery procedure. The procedure is straightforward to me,
> however its implementation in deal.ii not (yet). There are some functions like
> the VectorTools::project() function, which is supposed to give a L2-projections.
> However, I am not quite sure how to use them. Therefore I tried to implement
> the  Zienkiewicz-Zhu recovery procedure myself. But there is a fundamental
> problem I could not solve so far. In elasticity problems for 2D or 3D we have
> a vector-valued solution. Therefore, in 2D for a triangulation with n Nodes we
> have 2n DoFs. If I want to add the projected scalar function (e.g. the mises
> stress) in form of a command like
>
> data_out.add_data_vector(mises_stress_recovery,"mises_recovery");
>
> to the output, there is the problem that the mises_stress_recovery vector only
> has a dimension of n, where as data_out wants a vector of dimension 2n. How to
> deal with that issue?

The ZZ procedure creates a finite element field that represents the stress (or
one of its invariants). You need to say what finite element you want to use
for this, and create a DoFHandler to represent it. You do this here:

> I tried to generate a new dof_handler via
> DoFHandler<dim> dof_handler_recovery(this->triangulation);
> FESystem<dim> fe_recovery(FE_Q<dim>(this->fe_degree),1);// here dim=1 since we
> only have on dof per node
> MappingQ<dim> mapping_recovery(this->fe_degree);
> AffineConstraints<double> constraints_recovery;
> dof_handler_recovery.distribute_dofs(fe_recovery);
>
> and then assemble a mass matrix and a rhs to solve for the
> mises_stress_recovery vector. However, when looping over the cells, I am not
> sure how to reference to the cells. Should I use the cells of the newly
> created DoFHandler via

> for(constauto&cell :dof_handler_recovery.active_cell_iterators())

> or via the DoFHandler of the original problem solution via
> for(constauto&cell :dof_handler.active_cell_iterators()).

You have to have both. When you integrate over shape functions of the
recovered stress, you need to loop over the cells of dof_handler_recovery.
Where you need to evaluate the original solution, you need a cell of
dof_handler. There are functions that allow you to convert an iterator for a
cell to one DoFHandler to the corresponding iterator for a cell to another
DoFHandler.

> Do I really have to create a new DoFHandler? If not, how would one assemble
> the mass matrix, since for that we need the an local_dof_indices object like
> std::vector<types::global_dof_index> that points only to every second dof
> index or so.
> I would be grateful for some tips or directions.

If you don't want to use a different DoFHandler, then you need to add another
scalar field to the finite element system you are already using for the
original DoFHandler. That of course has consequences for the size of the
linear systems you build, so that may or may not be simpler.

[Konrad Schneider]

Wolfgang,
thanks again for your valuable advise. Always a pleasure to receive support from the boss ;-)

 

> In the scope of this loop the shape function gradients at the quadrature
> points are used to compute the desired output variable by a special
> calculation rule. However, these are not the vertices. If the gradients are
> somehow again computed at the vertices, which might be an internal procedure
> of the DataPostprocessor-functionality of dealii, how would this procedure
> know what the special calculation rule of the desired output would look like?

Can you clarify this question? What are the "special rules"?

Here I used the idiom "special rule" as a synonym for any post processing related computation such as computing the  von Mises stress, computing the maximum eigenstresses, evaluating the stress triaxiality etc. Basically any postprocessing quantity one can think.  However, since you explained, that the loop I was asking about, loops over evaluations points and not quadrature points, my question is answered.

 

You have to have both. When you integrate over shape functions of the
recovered stress, you need to loop over the cells of dof_handler_recovery.
Where you need to evaluate the original solution, you need a cell of
dof_handler. There are functions that allow you to convert an iterator for a
cell to one DoFHandler to the corresponding iterator for a cell to another
DoFHandler.

Following your encouragement of using two individual dof_handlers, I managed to do what I wanted following step 31. For completeness my code looks like the following:

// mises stress via ZZ - stress recovery

DoFHandler<dim> dof_handler_recovery(this->triangulation);

Vector<double> mises_stress_recovery(dof_handler_recovery.n_dofs());

{

// create new dof_handler

FE_Q<dim> fe_recovery(this->fe_degree);

MappingQ<dim> mapping_recovery(this->fe_degree);

AffineConstraints<double> constraints_recovery;

dof_handler_recovery.distribute_dofs(fe_recovery);

constraints_recovery.clear();

// setup system

DynamicSparsityPattern dsp(dof_handler_recovery.n_dofs());

DoFTools::make_sparsity_pattern(dof_handler_recovery,

dsp,

constraints_recovery,

/*keep_constrained_dofs = */ true);

SparsityPattern sparsity_pattern_recovery;

sparsity_pattern_recovery.copy_from(dsp);

SparseMatrix<double> mass_matrix;

mass_matrix.reinit(sparsity_pattern_recovery);

Vector<double> rhs_recovery;

rhs_recovery.reinit(dof_handler_recovery.n_dofs());

// assemble system

QGauss<dim> quadrature_formula(this->fe.degree + 1);

const unsigned int n_q_points = quadrature_formula.size();

FEValues<dim> fe_values_recovery(mapping_recovery,

fe_recovery,

quadrature_formula,

update_values | update_quadrature_points | update_JxW_values );

FEValues<dim> fe_values(this->mapping,

this->fe,

quadrature_formula,

update_gradients ); // since equal quadrature rule is used for both fe_values -> we only need to update the gradients

const unsigned int n_dof_per_cell = fe_recovery.n_dofs_per_cell();

FullMatrix<double> cell_matrix(n_dof_per_cell, n_dof_per_cell); // element mass matrix

Vector<double> cell_rhs(n_dof_per_cell); // element rhs-vector

std::vector<types::global_dof_index> local_dof_indices(n_dof_per_cell); // coincidence matrix local dofs -> global dof

const SymmetricTensor<4, 3> C_0=get_stress_strain_tensor<3>(this->E,this->nu);

auto cell = dof_handler_recovery.begin_active();

const auto endc = dof_handler_recovery.end();

auto stress_cell = dof_handler.begin_active();

for(;cell!= endc;++cell, ++stress_cell){

fe_values_recovery.reinit(cell);

fe_values.reinit(stress_cell);

cell_matrix = 0;

cell_rhs = 0;

const unsigned int n_dof_per_local_cell = cell->get_fe().dofs_per_cell;

std::vector<std::vector<Tensor<1, dim>>> u_gradients(n_q_points, std::vector<Tensor<1, dim>>(dim));

fe_values.get_function_gradients(this->u_glob, u_gradients);

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

for (unsigned int i=0; i<n_dof_per_local_cell; ++i){

for (unsigned int j=0; j<n_dof_per_local_cell; ++j){

cell_matrix(i,j) +=(fe_values_recovery.shape_value(i,q) *

fe_values_recovery.shape_value(j,q) *

fe_values_recovery.JxW(q));

}

// right-hand-side

Tensor<2, dim> grad_u;

for (unsigned int d = 0; d < dim; ++d) {

grad_u[d] = u_gradients[q][d]; // vector valued assignment

}

const SymmetricTensor<2, dim> strain_dim = symmetrize(grad_u);

SymmetricTensor<2, 3> strain;

if (dim<3){// plain strain case

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

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

strain[k][l] = strain_dim[k][l];

}

const SymmetricTensor<2, 3> stress= C_0 * strain;

const double mises_stress_qi = std::sqrt(3./2.)*deviator(stress).norm();

cell_rhs(i) += ( fe_values_recovery.shape_value(i,q) *

mises_stress_qi *

fe_values_recovery.JxW(q));

}

}

local_dof_indices.resize(cell->get_fe().dofs_per_cell);

cell->get_dof_indices(local_dof_indices);

constraints_recovery.distribute_local_to_global(cell_matrix,cell_rhs,local_dof_indices,mass_matrix,rhs_recovery);

}

// Solve the system

SparseDirectUMFPACK A_direct;

A_direct.initialize(mass_matrix);

A_direct.vmult(mises_stress_recovery, rhs_recovery);

}

data_out.add_data_vector(dof_handler_recovery,mises_stress_recovery, "mises_recovery");

 and the output is smooth as a baby's bottom ;-)

I also tried out the transfer properties of dof_handler cells you pointed out. Following the DofHandler-Question of the FAQ the cell-loop then has to start like this:

for(const auto & tria_cell : this->triangulation.active_cell_iterators()){

typename DoFHandler<dim>::active_cell_iterator cell (&this->triangulation, tria_cell->level(), tria_cell->index(), &dof_handler_recovery);

typename DoFHandler<dim>::active_cell_iterator stress_cell (&this->triangulation, tria_cell->level(), tria_cell->index(), &this->dof_handler);

However, in contrast to the FAQ answer one has to use the "typename" keyword.

Many thanks again and best regards
Konrad