Re: [deal.II] Neumann vector conditions

[luca.heltai]

> On 17 Feb 2017, at 14:46, Franco Milicchio  wrote:
> 
> I agree I have a translation problem here. In practice, I am integrating a 
> weak problem as this (in Fenics' UFL terms, but easy to recognize):
> 
> 
> inner( sigma(u), eps(v) ) * dx = inner(f, v) * dx + inner(f1, v) * ds(1)


I’m guessing that 

f1 = sigma(u) n 

is the Neumann data on the boundary. Then if you have f1 as a dim-dimensional 
vector field, it already contains n. 

the translation from ufl to deal.II is easy:


ds(1) =>  fe_face_values.JxW(f_q_point); 
v => fe_face_values.shape_value(i, f_q_point)
f1 =>  rhs_face_values[f_q_point]

since fe_face_values.shape_value(i, f_q_point) refers to the component_i of v, 
then

inner(f1, v) *ds(1) = sum over q of 

fe_face_values.shape_value(i, f_q_point) * 
rhs_face_values[f_q_point][component_i]) * fe_face_values.JxW(f_q_point); 


> where f is a bulk force (zero in my code), f1 is a (distributed) force 
> applied to the top boundary. The dx and ds terms refer to the integral on the 
> domain and boundary; u is the displacement, and v is the test function. Sigma 
> and eps are functions that compute, obviously, stress and strain tensors.
> 
> As I said, Fenics hid everything as you see above, but I am very eager to let 
> it go as soon as possible.
> 
>  
> 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); 
> 
> 
> My question here is this: why am I not using the normal to the face?

You are not using it, because you are imposing

sigma(u)n = f

on the boundary (where f is a dim-dimensional vector field on the boundary). 
Even in your FENICS implementation you are not using the normal… 

inner(f1,v)*ds(1) 

does not contain the normal anywhere, only the integral of f_i v_i on the faces 
on the boundary (I guess with id 1).


By the way this is the same for the Poisson problem. If you impose 

u,i n_i = f

then your integral on the boundary becomes

f*v*ds(1)












-- 
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.

Re: [deal.II] Neumann vector conditions

[Franco Milicchio]

On Friday, February 17, 2017 at 1:48:54 PM UTC+1, Luca Heltai wrote:
>
> Dear Franco, 
>
> I think there is a problem in your formulation… 
>
> You are integrating 
>
> u_i n_i f_i 
>
>
I agree I have a translation problem here. In practice, I am integrating a 
weak problem as this (in Fenics' UFL terms, but easy to recognize):


inner( sigma(u), eps(v) ) * dx = inner(f, v) * dx + inner(f1, v) * ds(1)


where f is a bulk force (zero in my code), f1 is a (distributed) force 
applied to the top boundary. The dx and ds terms refer to the integral on 
the domain and boundary; u is the displacement, and v is the test function. 
Sigma and eps are functions that compute, obviously, stress and strain 
tensors.

As I said, Fenics hid everything as you see above, but I am very eager to 
let it go as soon as possible.

 

> 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); 
>
>
My question here is this: why am I not using the normal to the face?

 

> 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. 
>


Thanks for your help, Luca!
 Franco

-- 
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.

Re: [deal.II] Neumann vector conditions

[luca.heltai]

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  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::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.