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
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/edbc2f13-176f-47f0-85e5-e3a0bd625ab0n%40googlegroups.com.
