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.

Reply via email to