there was a small bug in my code I posted above. the following code is
running fine.
Konrad Schneider schrieb am Freitag, 17. März 2023 um 17:03:18 UTC+1:
> 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/c0377cea-0a24-467f-bca9-c93be89dd8ben%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 DataPostprocessor<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(1);
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;
}
