Apply Neumann Boundary condition

[benhour....@gmail.com]

Dear All,

I want to apply non-zero neumann boundary condition on the curved domain of half of the circle. The way that I define my geometry, the boundary id and apply the Neumann B.C is written as follow:

Defining boundary id for different sections of domain (half_hyper_ball):

        const double tol_boundary = 1e-12

for (typename Triangulation<dim> ::active_cell_iterator

cell=triangulation.begin_active();

cell!=triangulation.end(); ++cell)

for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)

if (cell->face(face)->at_boundary() == true)

 {

   const Point<dim> face_center = cell->face(face)->center();

   if (std::abs(std::sqrt(face_center[0]*face_center[0]+face_center[1]*face_center[1])

                  - radius) < tol_boundary)

              cell->face(face)->set_boundary_id (2); // faces on the outer curved edge of the domain...

   else

    cell->face(face)->set_boundary_id (1);

// Implying the mechanical boundary value to the curved surface of the domain that appears on the right hand side 

    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)

      if (cell->face(face)->at_boundary() == true

          && cell->face(face)->boundary_id() == 2)

       {

        scratch.fe_face_values_ref.reinit(cell, face);

        for (unsigned int f_q_point=0; f_q_point<n_q_points_f; ++f_q_point)

     {

      const double eta_orderM = local_quadrature_point_history[f_q_point].get_eta_orderM();

      for (unsigned int i=0; i<dofs_per_cell; ++i)

       {

        const unsigned int i_group =fe.system_to_base_index(i).first.first;

        if (i_group == u_dof)

     {

      const Tensor<1, dim> neumann_value = eta_orderM * scratch.fe_face_values_ref.normal_vector(f_q_point);

      const unsigned int component_i = fe.system_to_component_index(i).first;

      const double JxW = scratch.fe_face_values_ref.JxW(f_q_point);

      data.cell_rhs(i) += neumann_value[component_i] * scratch.fe_face_values_ref.shape_value(i, f_q_point) * JxW;

     }

       }

     }

       }

It would be very kind of you if you let me know whether this definition for the Neumann boundary condition on the curved side of my domain is true or not.

Bests,

Benhour

[Jean-Paul Pelteret]

Dear Benhour,

As I explained to someone else the day before yesterday, the point of this forum is not to provide a code review service without you having done due diligence on your side. Does this not work as expected? If this is the case, then how is the problem manifesting itself? Unless you have a specific issue pertaining to the code that you have posted then I'm afraid that there's no discussion to be had here.

Regards,

Jean-Paul

[benhour....@gmail.com]

Dear Jean,

I know that this website is not the code reviewer. I want to apply Neumann boundary condition on the curved side of my domain. It should be note that the domain is symmetry. I have considered half of a circle. It would be very kind of you if you let me know whether it is applied correctly because I think something is wrong when I run my code and my domain completely distorted. I am doubtful about how I applied the boundary condition. For the horizontal axis, I have considered the zero displacement along y-axis because of symmetry. Is it correct? Should I fix the center of the circle as well? How Can I do that? In other words, How can I fix one point (vertex) in 2-D? Looking forward to hearing from you. 

[Xin-Bo Qi]

You should test your code step by step, figuring out exactly where the error occurs. And then propose a technical and accurate question here.

Time is expensive, and the developers offer their kind help to us, then we should respect their time and help.

在 2016年12月2日星期五 UTC+8上午7:25:32，benhour....@gmail.com写道：

[Jean-Paul Pelteret]

Dear Benhour,

Thank you for taking my point heart and providing more specifics about your actual problem.

For the horizontal axis, I have considered the zero displacement along y-axis because of symmetry. Is it correct? Should I fix the center of the circle as well?

If you don't provide at least one Dirichlet constraint for each of the solution components then your problem is indeterminate. So if you only constrain the x-DoFs along the y-axis and have no other Dirichlet constraints then this could be the source of some troubles. 

 

How Can I do that? In other words, How can I fix one point (vertex) in 2-D? 

You have to manually add such a point constraint to the constraint matrix or map of boundary values. If your centre point is coincident with a vertex of the triangulation then you can achieve this with the help of the cell->vertex_dof_index() function. There have been plenty of posts in the past on how this function works, so I need not explain it again in detail.

I should also note that for highly nonlinear problems, if your load increment is too large between time steps then this could also be problematic (unless you're using a load-following algorithm like the arc-length method). It may be of use to reduce the magnitude of the applied traction while you're debugging this.

Regards,

Jean-Paul

[benhour....@gmail.com]

Dear Jean,

Thanks very much for your help and support. I do really appreciate your time and know you are really busy. The fact is that my code ran just for the first time step and can not solve for the other time step. when I refine the global mesh twice, this problem happened, however, with increasing the number of global refine mesh, the whole geometry is distorted. In fact, it is a thermo-elastic problem which I should use two different Neumann B.C for the curved side of my domain. I just want to know whether applying this type was true or not. In addition, could you please help me more how I should use load-following algorithm in dealii? It should be noted that I have multiplied the Neumann magnitude by time-ramp (current_time/time_step). Looking forward to hearing from you.

[Wolfgang Bangerth]

On 12/02/2016 09:30 AM, benhour....@gmail.com wrote:
> Thanks very much for your help and support. I do really appreciate your time
> and know you are really busy. The fact is that my code ran just for the first
> time step and can not solve for the other time step. when I refine the global
> mesh twice, this problem happened, however, with increasing the number of
> global refine mesh, the whole geometry is distorted. In fact, it is a
> thermo-elastic problem which I should use two different Neumann B.C for the
> curved side of my domain. I just want to know whether applying this type was
> true or not. In addition, could you please help me more how I should use
> load-following algorithm in dealii? It should be noted that I have multiplied
> the Neumann magnitude by time-ramp (current_time/time_step). Looking forward
> to hearing from you.

Benhour -- it is practically impossible to help you with this question. You
mention the following topics:
* It is a time dependent problem
* There is adaptivity
* There seems to be mesh movement
* The geometry seems to be distorted
* It is a thermo-elastic problem
* There seem to be different boundary conditions
* There seem to be curved boundaries
* The boundary conditions are modified in some time dependent way
* There is mention of a "load-following" algorithm

These are 9 concepts, in 9 lines of text, all without any kind of elaboration.
In other words, we have no real idea what you are doing, and equally
important, we have no idea what is actually wrong. When debugging problems
like yours, you need to learn to reduce things to the minimal case that shows
the problem, and to look for the *first sign* where something seems to be
going wrong. In your case, if you know that the mesh is distorted, it is no
surprise that you cannot solve. So *why* is the mesh distorted? Have you
inspected the mesh before and after you move it? Have you inspected the
solution of the first time step that is presumably used to move the mesh?

These are the questions you need to ask yourself.

Best
W.

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

[Jean-Paul Pelteret]

Dear Benhour,

 

These are 9 concepts, in 9 lines of text, all without any kind of elaboration.
In other words, we have no real idea what you are doing, and equally
important, we have no idea what is actually wrong. When debugging problems
like yours, you need to learn to reduce things to the minimal case that shows
the problem, and to look for the *first sign* where something seems to be
going wrong.

Prof Bangerth has succinctly stated why I'm having such difficulties in helping you. 

I, like him, would encourage you to take your code and strip out as many of the complexities as possible, and produce some test cases that convince you that each part of your code is working as expected. This is exactly how I build up any of my codes that tackle complex problems. For this you often need to test your implementation of each "feature" against a really simple problem that you can hand calculate a comparison result for, or otherwise some benchmark in the literature.

Thats really the best advice we can give you in this instance. So, for example, coming back to your question of the implementation of the traction boundary condition: you could easily test this by considering a distributed load on a bar (inducing uni-axial tension) or perhaps the result shown in step-44 or one of the more simple code-gallery examples.

Regards,

Jean-Paul

[benhour....@gmail.com]

Dear Jean,

Thanks very much for your help and support. Most of my problem is related to applying axisymmetry boundary condition for a curved domain. I attached my code for your consideration. According to the fact that this is my first project with dealii, Please accept my lack of knowledge in this case. This is a themo-elastic problem for a curved domain which for getting result, I have fixed the straight line of my domain. However, I really do not know why the Neumann boundary condition which is brought in the assemble-system-rhs does not work properly, instead of putting force on the curved domain, it seems that it is a tensile force. It should be noted that I have changed the sign of my equation for this type of boundary condition, however, it  still has not worked properly. Your noteworthy comments on my code can give me incredible experiences on working with dealii. Looking forward to hearing from you.

Bests,

Benhour

[Wolfgang Bangerth]

On 12/06/2016 08:26 AM, benhour....@gmail.com wrote:
> Thanks very much for your help and support. Most of my problem is related to
> applying axisymmetry boundary condition for a curved domain. I attached my
> code for your consideration. According to the fact that this is my first
> project with dealii, Please accept my lack of knowledge in this case. This is
> a themo-elastic problem for a curved domain which for getting result, I have
> fixed the straight line of my domain. However, I really do not know why the
> Neumann boundary condition which is brought in the assemble-system-rhs does
> not work properly, instead of putting force on the curved domain, it seems
> that it is a tensile force. It should be noted that I have changed the sign of
> my equation for this type of boundary condition, however, it still has not
> worked properly. Your noteworthy comments on my code can give me incredible
> experiences on working with dealii. Looking forward to hearing from you.

Benhour, please do read my and J-P's emails again and think how they may apply
to the situation you describe. Your problem is still using (i) curved domains,
(ii) axisymmetric coordinate systems, (iii) couples multiple equations, (iv)
tries to solve a complicated first before testing a simple case, and numerous
other things.

We are, by and large, very happy to answer specific questions on this forum.
But you can't expect that people here are willing to take the 2238 lines of
code you sent and debug the problem for you.

Please reconsider what J-P and I said in previous emails.