On Thursday, February 16, 2017 at 7:49:52 PM UTC+1, Jean-Paul Pelteret
wrote:
>
> 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>
> .
>
Thanks for your pointer, Jean-Paul, I appreciate it. I am coming from
Fenics, where all these details were completely hidden and almost
inaccessible: I just used weak forms there. It was easier but, ultimately,
not powerful enough. So this is new to me, I am sorry if I'm asking
something trivial.
I have tried to implement it, but I think I have a wrong result. This is
the code I've written:
// Neumann (boundary forces on top boundary == id 2)
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))
{
// Get face FE values for this face and cell
fe_face_values.reinit(cell, face_number);
// Neumann force vector on the face
right_face_side.vector_value_list(fe_face_values.
get_quadrature_points(), rhs_face_values);
for(unsigned int f_q_point = 0; f_q_point <
n_face_q_points; ++f_q_point)
{
// Add contributions to all DOFs
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const unsigned int component_i = fe.
system_to_component_index(i).first;
cell_rhs(i) += ( fe_face_values.shape_value(i,
f_q_point) *
rhs_face_values[f_q_point][
component_i] *
fe_face_values.normal_vector
(f_q_point)[component_i]
) * fe_face_values.JxW(f_q_point);
}
} // Integral for linear form
}
} // Boundary forces
My point is to integrate a distributed force f on the top boundary, so as
far as I understand this (from the tutorial 44), I have to do cycle over
all boundary faces, get values (normals, shape functions, jacobians) and
their values on the quadrature points, and finally add them to the right
hand side.
Where is the error?
Thank you very much!
Franco
PS. I am still struggling on some terms in deal.II, as for instance the
role of "fe.system_to_component_index(i).first", but I need to write code
to understand thoroughly how to switch my codebase from Fenics.
--
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.
