One more thing:
you should loop over q point before loop over dof.

Léonhard YU <[email protected]> 于2024年6月8日周六 18:14写道:

> Dear all,
>  I am working on impose constant pressure load on a plate's upper surface.
> The code is derived from step-18 and attached below. But I found that the
> diplacement is not right, and it varies when the number of face elements
> changes. I am thinking about something wrong with the loading but I do not
> have any idea. Thanks a lot if you can help me.
> The constant pressure is 1e-2 and the boundary_id for loading is
> bc_uppersurface. The codes attached below:
>
> for (const auto &cell : dof_handler.active_cell_iterators())
>     if (cell->is_locally_owned())
>     {
>         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); //
>             }
>
>     // pressure load
>         for (const auto &face : cell->face_iterators())
>         {
>             if (face->at_boundary() && (face->boundary_id() ==
> bc_uppersurface))
>             {
>             fe_face_values.reinit(cell, face);
>             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 face_q_point = 0; face_q_point <
> face_n_q_point; ++face_q_point)
>                 {
>                 Tensor<1,dim> dir;
>                 dir = fe_face_values.normal_vector(face_q_point);
>                 Point<dim> point_dof =
> fe_face_values.quadrature_point(face_q_point);
>                 const double magnitude = 1e-2;
>                 const Tensor<1, dim> traction = -1.0 * magnitude * dir;
>                 cell_rhs(i) +=
>                             fe_face_values.shape_value(i, face_q_point) *
> traction[component_i] *
>                         fe_face_values.JxW(face_q_point);
>                 }
>             }
>             }
>         }
>     }
>
> Best Regards
>
> --
> 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/d9ed8d3d-aa19-4367-9505-ca9d9656d45cn%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/d9ed8d3d-aa19-4367-9505-ca9d9656d45cn%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>

-- 
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/CADr3Od%2Bi2N7Ls8_J7S-aWWfcQ-D5XP8iLky4bS40yXegcm1xYA%40mail.gmail.com.

Reply via email to