Dear Franco,
Super quick answer: Step-44 demonstrates how to implement the Neumann BC
for elasticity
<https://www.dealii.org/developer/doxygen/deal.II/step_44.html#Solidassemble_system_rhs>
.
Best,
Jean-Paul
On Thursday, February 16, 2017 at 6:22:06 PM UTC+1, Franco Milicchio wrote:
>
> Dear all,
>
> I am implementing a very simple linear elastic problem, and I've succeeded
> in using Dirichlet conditions. However, reading the tutorial I am unable to
> implement a Neumann condition. In my simple example I have something like
> this:
>
> ^^^^^^^^^^ f
> +----------+
> | |
> | |
> | |
> | |
> +----------+
> //////////// clamped
>
> I don't understand how I can implement the integral on the top boundary,
> at the moment I've written this code (which does not compile):
>
> // Force (neumann bc), no force for now
>
> ConstantFunction<2> right_hand_side(std::vector<double>({ 0.0, 0.0
> }));
>
>
>
> // Force on boundary (neumann bc)
>
> ConstantFunction<2> right_face_side(std::vector<double>({ 0.0, 111.0
> }));
>
>
> // For each cell
>
> for (auto cell : dof_handler.active_cell_iterators())
>
> {
>
> // Zero local matrix and vector
>
> cell_matrix = 0;
>
> cell_rhs = 0;
>
>
>
> // Compute FE values for this cell
>
> fe_values.reinit(cell);
>
>
>
> // Compute lambda and mu (in case we can use custom functions)
>
> lambda.value_list(fe_values.get_quadrature_points(),
> lambda_values);
>
> mu.value_list(fe_values.get_quadrature_points(), mu_values);
>
>
> // Neumann force vector
>
> right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
> rhs_values);
>
>
> // ...assemble lhs matrix and rhs vector
>
>
> // Neumann (boundary forces, 2 is top boundary)
>
> for (unsigned int face_number = 0; face_number < GeometryInfo
> <dim>::faces_per_cell; ++face_number)
>
> {
>
> if (cell->face(face_number)->at_boundary() && (cell->face
> (face_number)->boundary_id() == 2))
>
> {
>
> const unsigned int n_face_q_points =
> face_quadrature_formula.size();
>
>
> fe_face_values.reinit(cell, face_number);
>
>
>
> auto qps = fe_face_values.get_quadrature_points();
>
>
> for (unsigned int q_point = 0; q_point < n_face_q_points;
> ++q_point)
>
> {
>
> const double neumann_value =
> right_face_side.value(fe_face_values.quadrature_point(q_point));
>
>
>
> for (unsigned int i = 0; i < dofs_per_cell; ++i)
>
> {
>
> cell_rhs(i) += (neumann_value *
>
> fe_face_values.shape_value(i,
> q_point) *
>
> fe_face_values.JxW(q_point));
>
> }
>
> }
>
> }
>
> } // Boundary forces
>
> I have also tried to use the normal vector, but weirdly enough, I don't
> get a tension, but a shear with the following code:
>
> const double neumann_value = (right_face_side.value
> (qps[q_point], 0) * fe_face_values.normal_vector(q_point)[0])
>
> + (right_face_side.value
> (qps[q_point], 1) * fe_face_values.normal_vector(q_point)[1]);
>
> How can I implement this?
>
> Thanks!
>
--
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.
