Applying non-homogeneous Neumann BC in step-18

[Mohammad Amir Kiani Fordoei]

Dear Deal.II group

I am trying to add boundary force (b) to step-18 as below

As I understood, I just need to specify a boundary_id on triangulation and apply traction force on that with code like this:
{ - loop over cells
-loop over faces per cell && controlling being at_boundary
-loop over quadrature points of face
-loop over dof per cell && controlling intended boundary_id_NBC
-finally adding this term to rhs as
cell_rhs (i)+=

(fe_face_values.shape_value (i,q_point)

*traction_component

*fe_face_values.JxW(q_point));}

But, I didn't see noticble variation of displacement or norm of stress in my solution, even for large traction force.
is it necessary to multiply traction_component to "* present_timestep * velocity " (like Dirichlet BC displacement)? (I tried this too, but similar to previous one no significant effect of load was observed.)

instead of procedure described in step-18, can I use functions " create_right_hand_side" and "create_boundary_right_hand_side" of VectorTools for applying body and traction force, respectively or they just belongs to Laplace equation?
Thanks for considering this!

Amir

[Wolfgang Bangerth]

On 7/11/23 05:32, Mohammad Amir Kiani Fordoei wrote:
> Screenshot from 2023-07-11 14-01-30.png
> As I understood, I just need to specify a boundary_id on triangulation and
> apply traction force on that with code like this:
> { - loop over cells
> -loop over faces per cell && controlling being at_boundary
> -loop over quadrature points of face
> -loop over dof per cell && controlling intended boundary_id_NBC
> -finally adding this term to rhs as
> cell_rhs (i)+=
> (fe_face_values.shape_value (i,q_point)
>
> *traction_component
>
> *fe_face_values.JxW(q_point));}
> But, I didn't see noticble variation of displacement or norm of stress in my
> solution, even for large traction force.
> is it necessary to multiply traction_component to "* present_timestep *
> velocity " (like Dirichlet BC displacement)? (I tried this too, but similar to
> previous one no significant effect of load was observed.)

Is this a boundary where you are still *also* applying a prescribed
displacement? You can only do one or the other.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Mohammad Amir Kiani Fordoei]

Dear Wolfgang
Sorry for late response. Due to internet censorship in Iran, I didn't have access to google group for few days.
I consider a specific region on triangulation for Neumann BC (a force vector in 3D space).
There must be some mistakes in my program, but I couldn't figure out.

I tried to follow your instruction on:
https://groups.google.com/g/dealii/c/_4eeywhhxk4/m/myLiTWCqBgAJ .

So, I used this:
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int q_point = 0; q_point < n_face_q_points ; ++q_point)
{
const Tensor<1,dim> n = fe_face_values.normal_vector(q_point);
cell_rhs (i)+=( (n*traction_component*fe_face_values.shape_value(i, q_point))*
fe_face_values.JxW(q_point));
}//for q_point
}//for dof_per_cell

But in computation of cell rhs i get this error:

/home/amir/eclipse-workspace/mech1/mech.cc:983:22: error: no match for ‘operator+=’ (operand types are ‘double’ and ‘dealii::Tensor<1, 3>’)
  983 |          cell_rhs (i)+=( (n*traction_component*fe_face_values.shape_value(i, q_point))*
      |          ~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  984 |                   fe_face_values.JxW(q_point));

As error says, the normal vector with 3 component cannot be multiplied by double values in other terms.
Additionally,  my traction component (for example in x-direction) is preserved for dofs of other directions.
I know these mistakes, but I  cannot figure out proper solution for them.
I have read steps 7, 18, 20, 21 and also reviewed lectures 21.~21.65. But I couldn't find similar codes.

Thanks for your consideration.

Best regards

Amir

[Wolfgang Bangerth]

On 7/15/23 11:33, Mohammad Amir Kiani Fordoei wrote:
>
> const Tensor<1,dim> n = fe_face_values.normal_vector(q_point);
> cell_rhs (i)+=( (n*traction_component*fe_face_values.shape_value(i, q_point))*
> fe_face_values.JxW(q_point));
> }//for q_point
> }//for dof_per_cell
>
> But in computation of cell rhs i get this error:
>
> /home/amir/eclipse-workspace/mech1/mech.cc:983:22: error: no match for
> ‘operator+=’ (operand types are ‘double’ and ‘dealii::Tensor<1, 3>’)
>   983 |          cell_rhs (i)+=(
> (n*traction_component*fe_face_values.shape_value(i, q_point))*
>       |
>  ~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>   984 |                   fe_face_values.JxW(q_point));
>
> As error says, the normal vector with 3 component cannot be multiplied by
> double values in other terms.

No, that's not actually what the error message says :-) It doesn't say
anything about multiplication. If you read it carefully, it says that you
cannot assign a right hand side of type Tensor<1,3> to a left hand side of
type double. That makes sense.

The types you have are because your right hand side is the product of
n -- a tensor
traction_component -- not sure, but probably a scalar
fe_face_values.shape_value(i, q_point)) -- a scalar, see the function's
documentation
fe_face_values.JxW(q_point) -- a scalar

What you *want* is a vector-valued shape function. You get that by
"extracting" this from the fe_values object in the same way as is discussed in
the vector-valued documentation module. I think something like this should work:

FEValuesExtractors::Vector displacement(0);
cell_rhs (i)+=(
(n*traction_component*fe_face_values[displacement].shape_value(i, q_point))*
fe_face_values.JxW(q_point));