I tried to use code from step-22 and from the sources you recommended.
The problem is that the solution is zero over the body except the boundary
where
Dirichlet conditions were set. It is just one iteration. I don't know why
they are not in use.
Briefly, my code consists of the following steps which are perfomed in turn
(from 1 to 5)
1) Setup variables and constraints
dof_handler.distribute_dofs (fe);
{
constraints.clear ();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(dof_handler,
1,
boundary_u_handler[1],
constraints);
VectorTools::interpolate_boundary_values(dof_handler,
3,
boundary_u_handler[3],
constraints);
}
constraints.close ();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from (dsp);
system_matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
2) Assembling FEM matrix
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.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);
// Loop over all cells:
for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end();
++cell) {
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit(cell);
for (unsigned int i = 0; i < dofs_per_cell; ++i) {
for (unsigned int j = 0; j < dofs_per_cell; ++j) {
for (unsigned int q_point = 0; q_point < n_q_points;
++q_point) {
const SymmetricTensor<2, dim>
eps_phi_i = get_strain(fe_values, i, q_point),
eps_phi_j = get_strain(fe_values, j, q_point);
cell_matrix(i, j) +=
(eps_phi_i * stress_strain_tensor * eps_phi_j
*
fe_values.JxW(q_point));
}
}
}
cell->get_dof_indices (local_dof_indices);
constraints.
distribute_local_to_global(cell_matrix, cell_rhs,
local_dof_indices,
system_matrix, system_rhs);
}
3) Assembling RHS vector
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
update_values |
update_normal_vectors |
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
// right hand side vector for the cell
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
std::vector<Vector<double> > rhs_values (n_q_points,
Vector<double>(dim));
// reinit right hand side
system_rhs.reinit (dof_handler.n_dofs());
// loop over all cells:
for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end();
++cell) {
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
// Assembling the right hand side boundary terms
for (unsigned int face_number=0;
face_number<GeometryInfo<dim>::faces_per_cell; ++face_number) {
Tensor<1, dim> neumann_vector_value;
if (cell->face(face_number)->at_boundary()) {
int id = cell->face(face_number)->boundary_id();
if (boundary_force.find(id) != boundary_force.end()) {
fe_face_values.reinit(cell, face_number);
// quadrature
for (unsigned int i = 0; i < dofs_per_cell; ++i) {
const unsigned int component_i =
fe.system_to_component_index(i).first;
for (unsigned int q = 0; q < n_face_q_points; ++q) {
boundary_force_handler.force_value(
fe_face_values.quadrature_point(q),
neumann_vector_value,
fe_face_values.normal_vector(q),
id);
cell_rhs(i) +=
neumann_vector_value[component_i] *
fe_face_values.shape_value(i, q)
*
fe_face_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);
}
4) Solve
SolverRichardson<> rich(solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
rich.solve (system_matrix, solution, system_rhs, preconditioner);
5)
constraints.distribute (image_solution);
The problem is that left boundary is fixed, up and down is under zero
force, and right boundary is displaced to the right.
The solution without 5) is zero, and with 5) it is zero everywhere outside
non-zero dirichlet boundary.
However the correct solution is linear x_displacement profile from zero on
the right boundary to non-zero value on the left.
I uploaded figure in the post above.
Regards
On Thursday, February 22, 2018 at 5:00:45 PM UTC+3, Denis Davydov wrote:
>
> Hi IIiya,
>
> On Thursday, February 22, 2018 at 11:55:33 AM UTC+1, Bryukhanov Ilya wrote:
>>
>> Denis, thanks a lot for the answer!
>>
>> I think that on the first step I can use usual
>> "interpolate_boundary_values + apply_boundary_values" functions,
>>
>
> the outline I gave above is agnostic to "first step"/"other steps", you
> would call VectorTools::interpolate_boundary_values() regardless
> stage to get constraints into the ConstraintMatrix object. If you do
> Newton-Raphson and solve for increments, at the first non-linear step you
> should have all constraints (hanging nodes + non-zero Dirichlet for a given
> step in quasi-static(?) loading), and for consequent steps only hanging
> nodes plus zero Dirichlet.
>
> Please read this module:
> https://www.dealii.org/developer/doxygen/deal.II/group__constraints.html
>
> This code-gallery may help in figuring out how to do this:
> https://www.dealii.org/developer/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html
>
>
>
>
>> that transform the global matrix: columns of the constrained nodes become
>> zero columns with diagonal non-zero entries.
>>
>> In a time loop I have to account boundary condition without modifing
>> global matrix. While looping over dirichlet constrained cells
>> I call the function
>> "hanging_node_constraints.distribute_local_to_global(cell_rhs,
>> local_dof_indices, system_rhs, cell_matrix)"
>> which modifies the system_rhs vector by adding columns from constrained
>> nodes, which are actually zero in global matrix.
>>
>
> Nope, if your "hanging_node_constraints" indeed contains only constraints
> from hanging nodes, you won't see any contribution to the RHS
> from it as constraints are homogeneous (ie. you wont' have u_{123} = 456
> from hanging nodes).
>
>
>>
>> Also I need to add components to rhs vector from the force boundary
>> conditions.
>>
>
> that's a different story, just do integration over faces and distribute as
> usual.
>
>
>>
>> However, I have to somehow account Dirichlet boundary conditions as in
>> previous step I didn't use any boundary values. The function
>>
>
> use ConstraintMatrix which has non-zero Dirichlet BC constraints.
>
>
>> "interpolate_boundary_values(dof_handler, boundary_id, DirichletHandler,
>> boundary_values)"
>> gives me "boundary_values" variables that maps node numbers to dirichle
>> boundary.
>>
>> But I don't know how should I use it. What should I do with that variable
>> to account boundary conditions together with "distribute_local_to_global"?
>>
>> Thanks a lot in advance. Sorry for my misunderstanding.
>>
>
> Regards,
> Denis.
>
--
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.
For more options, visit https://groups.google.com/d/optout.
