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 [email protected].
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;
}
