Dear Franco,
I think there is a problem in your formulation…
You are integrating
u_i n_i f_i
and while you can certainly come up with good reasons for this to work out, I’m
unsure about what you want to achieve. If your boundary conditions are
((grad(u)).n,v) = u_i,j n_j v_i
and you impose
(grad(u).n)_i = f_i
on the boundary, then you should have
cell_rhs(i) += ( fe_face_values.shape_value(i, f_q_point) *
rhs_face_values[f_q_point][component_i]) *
fe_face_values.JxW(f_q_point);
on the other hand, if you want to have, for example, for a scalar function p,
u_i,j n_j = p n_i
then you should have
cell_rhs(i) += ( fe_face_values.shape_value(i, f_q_point) *
rhs_scalar_function[f_q_point]
fe_face_values.normal_vector(f_q_point)[component_i]
) * fe_face_values.JxW(f_q_point);
What is the exact term you want to integrate?
L.
> On 17 Feb 2017, at 13:03, Franco Milicchio <franco.milicc...@gmail.com> wrote:
>
>
>
> 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.
>
>
> 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.
--
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.
