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 -- 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/c711eb18-b84b-4120-a46e-5bfc76773ec0n%40googlegroups.com.
#include <deal.II/base/quadrature_lib.h> #include <deal.II/base/function.h> #include <deal.II/base/tensor.h> #include <deal.II/lac/vector.h> #include <deal.II/lac/full_matrix.h> #include <deal.II/lac/sparse_matrix.h> #include <deal.II/lac/dynamic_sparsity_pattern.h> #include <deal.II/lac/solver_cg.h> #include <deal.II/lac/precondition.h> #include <deal.II/lac/affine_constraints.h> #include <deal.II/grid/tria.h> #include <deal.II/grid/grid_generator.h> #include <deal.II/grid/grid_refinement.h> #include <deal.II/grid/grid_out.h> #include <deal.II/dofs/dof_handler.h> #include <deal.II/dofs/dof_tools.h> #include <deal.II/fe/fe_values.h> #include <deal.II/numerics/vector_tools.h> #include <deal.II/numerics/matrix_tools.h> #include <deal.II/numerics/data_out.h> #include <deal.II/lac/sparse_direct.h> #include <deal.II/fe/fe_system.h> #include <deal.II/fe/fe_q.h> #include <fstream> #include <iostream> double TOL = 10E-9; namespace Step8{ using namespace dealii; template <int dim> class ElasticProblem{ public: ElasticProblem(); void run(); private: double E=10000.,nu=0.3; void setup_system(); void assemble_system(); void solve(); void output_results() const; Triangulation<dim> triangulation; DoFHandler<dim> dof_handler; FESystem<dim> fe; AffineConstraints<double> constraints; SparsityPattern sparsity_pattern; SparseMatrix<double> system_matrix; Vector<double> solution; Vector<double> system_rhs; }; template <int dim> void right_hand_side(const std::vector<Point<dim>> &points,std::vector<Tensor<1, dim>> & values){ AssertDimension(values.size(), points.size()); Assert(dim >= 2, ExcNotImplemented()); Point<dim> point_1, point_2; point_1(0) = 0.5; point_2(0) = -0.5; for (unsigned int point_n = 0; point_n < points.size(); ++point_n) { if (((points[point_n] - point_1).norm_square() < 0.2 * 0.2) || ((points[point_n] - point_2).norm_square() < 0.2 * 0.2)) values[point_n][0] = 1.0; else values[point_n][0] = 0.0; if (points[point_n].norm_square() < 0.2 * 0.2) values[point_n][1] = 1.0; else values[point_n][1] = 0.0; } } template <int dim> ElasticProblem<dim>::ElasticProblem() : dof_handler(triangulation) , fe(FE_Q<dim>(1), dim) {} template <int dim> void ElasticProblem<dim>::setup_system(){ dof_handler.distribute_dofs(fe); solution.reinit(dof_handler.n_dofs()); system_rhs.reinit(dof_handler.n_dofs()); constraints.clear(); for (const auto &cell : triangulation.active_cell_iterators()){ for (const auto &face : cell->face_iterators()){ if (face->at_boundary()){ if (std::fabs(face->center()[0] + 0.5) < TOL) face->set_boundary_id(100); else if (std::fabs(face->center()[1] + 0.5) < TOL) face->set_boundary_id(101); else if (std::fabs(face->center()[0] - 0.5) < TOL) face->set_boundary_id(103); else if (std::fabs(face->center()[1] - 0.5) < TOL) face->set_boundary_id(104); else face->set_boundary_id(10); } } } VectorTools::interpolate_boundary_values( dof_handler, // DOF handler object 100, // boundary id Functions::ZeroFunction<dim>(dim), // homogeneous boundary values constraints, // constraint object to be modified ComponentMask(std::vector<bool>{true, false})); VectorTools::interpolate_boundary_values( dof_handler, 101, Functions::ZeroFunction<dim>(dim), constraints, ComponentMask(std::vector<bool>{false, true})); VectorTools::interpolate_boundary_values( dof_handler, 103, Functions::ConstantFunction<dim>(std::vector<double>{0.5, 0.0}), // pull boundary values constraints, ComponentMask(std::vector<bool>{true, false})); constraints.close(); DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, /*keep_constrained_dofs = */ false); sparsity_pattern.copy_from(dsp); system_matrix.reinit(sparsity_pattern); } template <int dim> SymmetricTensor<4, dim> get_stress_strain_tensor(const double E, const double nu){ const double lambda = (E*nu)/((1.+nu)*(1.-2.*nu)); const double mu = E/(2.*(1.+nu)); SymmetricTensor<4, dim> C; for (unsigned int i = 0; i < dim; ++i) for (unsigned int j = 0; j < dim; ++j) for (unsigned int k = 0; k < dim; ++k) for (unsigned int l = 0; l < dim; ++l) C[i][j][k][l] = ( ((i == k) && (j == l) ? mu : 0.0) + ((i == l) && (j == k) ? mu : 0.0) + ((i == j) && (k == l) ? lambda : 0.0)); // ? returns value on the condition -> if condition is true return mu or lambda, else 0. return C; } template <int dim> inline SymmetricTensor<2, dim> get_B_tensor( const FEValues<dim> &fe_values,const unsigned int shape_func,const unsigned int q_point){ /* this is the B-Matrix or the B-vector in the integral K_ij^el = \int_{\Omega} \vec{B}_i \mat{C} \vec{B}_j dx which might be reffered to as strains since we have \int_{\Omega} \eta_{i,j} \sigma_{ij} dV = \int_{Omega} \eta_j t_j dA \int_{\Omega} \frac{1}{2}[\eta_{i,j} + \eta_{j,i}] \sigma_{ij} dV = \int_{Omega} \eta_j t_j dA \int_{\Omega} \varepsilon^{\eta}_{ab} \sigma_{ab} dV = \int_{Omega} \eta_j t_j dA with u_a = \varphi_{aj} u^F_j where u^F_j are the to be determined degrees of freedom and \eta_b = \varphi_{aj} \eta^F_j where \eta^F_j are the arbitraray values of the test functions we have \varepsilon_{ab} = \frac{1}{2}[u_{a,b}+u_{b,a}] = \frac{1}{2}[\varphi_{aj,b}+\varphi_{bj,a}] u^F_j \varepsilon_{ab} = B_{ajb} u^F_j and \varepsilon_{ab}^\eta = \frac{1}{1}[\eta_{a,b}+\eta_{b,a}] = = \frac{1}{2}[\varphi_{ai,b}+\varphi_{bi,a}] \eta^F_i \varepsilon_{ab}^\eta = B_{aib} \eta^F_i now we plugin these expressions into the integral with \sigma_{ab} = C_abmn and get \eta^F_i \int_{\Omega} B_{aib} C_{abmn} B_{mjn} dV u^F_j = \eta^F_i \int_{Omega} \varphi_{ai} t_a dA \eta^F_i K_ij u^F_j = \eta^F_i F_i which must hold for all \eta^F_i such that we get K_ij u^F_j = F_i This function implements the quantity B_{ikj}(x_q)=B_{ij} for a given shape function k associated with a dof and given x_q (quadrature point) such that B_{ij} is a second order Tensor w.r.t. i and j which can be used for the double contraction with the stiffness tensor C_{ijkl} */ SymmetricTensor<2, dim> B; // normal components for (unsigned int i = 0; i < dim; ++i) B[i][i] = fe_values.shape_grad_component(shape_func, q_point, i)[i]; // shear components for (unsigned int i = 0; i < dim; ++i) for (unsigned int j = i + 1; j < dim; ++j) B[i][j] = ( fe_values.shape_grad_component(shape_func, q_point, i)[j] + fe_values.shape_grad_component(shape_func, q_point, j)[i]) /2; return B; } template <int dim> void ElasticProblem<dim>::assemble_system(){ QGauss<dim> quadrature_formula(fe.degree + 1); FEValues<dim> fe_values(fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); const unsigned int n_q_points = quadrature_formula.size(); FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell); Vector<double> cell_rhs(dofs_per_cell); std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell); const auto C_0 = get_stress_strain_tensor<dim>(this->E,this->nu); for (const auto &cell : dof_handler.active_cell_iterators()){ cell_matrix = 0; cell_rhs = 0; fe_values.reinit(cell); const unsigned int n_dof_per_cell = cell->get_fe().dofs_per_cell; for (unsigned int q = 0; q < n_q_points; ++q){ for (unsigned int i=0; i<n_dof_per_cell; ++i){ for (unsigned int j=0; j<n_dof_per_cell; ++j){ const auto B_i = get_B_tensor<dim>(fe_values, i, q); const auto B_j = get_B_tensor<dim>(fe_values, j, q); cell_matrix(i,j) += ( B_i * C_0 * B_j )*fe_values.JxW(q); } } } cell->get_dof_indices(local_dof_indices); constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); } } template <int dim> void ElasticProblem<dim>::solve(){ SparseDirectUMFPACK A_direct; A_direct.initialize(this->system_matrix); A_direct.vmult(this->solution, this->system_rhs); constraints.distribute(solution); } // POST-PROCESSING template<unsigned int dim> class StressPostprocessor : public <dim>{ public: StressPostprocessor(double E,double nu); double E,nu; virtual void evaluate_vector_field( const DataPostprocessorInputs::Vector<dim> &inputs, std::vector<Vector<double>> &computed_quantities) const override; virtual std::vector<std::string> get_names() const override; virtual std::vector<DataComponentInterpretation::DataComponentInterpretation>get_data_component_interpretation() const override; virtual UpdateFlags get_needed_update_flags() const override; }; template <unsigned int dim> StressPostprocessor<dim>::StressPostprocessor(double E,double nu) : E(E) , nu(nu) {} template <unsigned int dim> std::vector<std::string> StressPostprocessor<dim>::get_names() const { std::vector<std::string> solution_names; solution_names.emplace_back("sig_11"); solution_names.emplace_back("sig_22"); solution_names.emplace_back("sig_33"); solution_names.emplace_back("sig_12"); solution_names.emplace_back("sig_13"); solution_names.emplace_back("sig_23"); solution_names.emplace_back("Mises"); return solution_names; } template <unsigned int dim> std::vector<DataComponentInterpretation::DataComponentInterpretation> StressPostprocessor<dim>::get_data_component_interpretation() const { std::vector<DataComponentInterpretation::DataComponentInterpretation> interpretation(7,DataComponentInterpretation::component_is_scalar); return interpretation; } template <unsigned int dim> UpdateFlags StressPostprocessor<dim>::get_needed_update_flags() const{ return update_gradients; } template <unsigned int dim> void StressPostprocessor<dim>::evaluate_vector_field(const DataPostprocessorInputs::Vector<dim> &inputs, std::vector<Vector<double>> &computed_quantities) const{ AssertDimension(computed_quantities.size(), inputs.solution_gradients.size()); const unsigned int n_quadrature_points = inputs.solution_values.size(); const auto C_0 = get_stress_strain_tensor<3>(this->E,this->nu); for (unsigned int q = 0; q < n_quadrature_points; ++q){ Tensor<2, dim> grad_u; for (unsigned int d = 0; d < dim; ++d) grad_u[d] = inputs.solution_gradients[q][d]; const SymmetricTensor<2, dim> strain_dim = symmetrize(grad_u); SymmetricTensor<2, 3> strain; if (dim<3){// plain strain 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; computed_quantities[q](0) = stress[0][0]; computed_quantities[q](1) = stress[1][1]; computed_quantities[q](2) = stress[2][2]; computed_quantities[q](3) = stress[0][1]; computed_quantities[q](4) = stress[0][2]; computed_quantities[q](5) = stress[1][2]; computed_quantities[q](6) = std::sqrt(3./2.)*deviator(stress).norm(); // mises stress } } template <int dim> void ElasticProblem<dim>::output_results() const{ DataOut<dim> data_out; data_out.attach_dof_handler (dof_handler); std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(dim, DataComponentInterpretation::component_is_part_of_vector); data_out.add_data_vector (solution, std::vector<std::string>(dim,"displacement"), DataOut<dim>::type_dof_data, data_component_interpretation); StressPostprocessor<dim> stress_post(this->E,this->nu); data_out.add_data_vector(solution, stress_post); data_out.build_patches (); std::ofstream output("solution.vtu"); data_out.write_vtu(output); } template <int dim> void ElasticProblem<dim>::run(){ const double half_edge_length = 0.5; const double inner_radius = 0.1; GridGenerator::hyper_cube_with_cylindrical_hole(triangulation, inner_radius, // inner_radius half_edge_length, // half the edge length of the square 0.5, // extension in z-dir 1, // repetion in z-dir false // colorize ); triangulation.refine_global(2); std::cout << " Number of active cells: " << triangulation.n_active_cells() << std::endl; setup_system(); std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; assemble_system(); solve(); output_results(); } } // namespace Step8 int main(){ try { Step8::ElasticProblem<2> elastic_problem_2d; elastic_problem_2d.run(); } catch (std::exception &exc) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Exception on processing: " << std::endl << exc.what() << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } return 0; }