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