Using p-refinement with high order elements

[Deepak Gupta]

Dear All,

I am trying to solve an elasticity problem and am dealing with p-refinement. With global refinement, my solution is as expected (correct). I am using an initial mesh of 4X2 elements. If I keep p=1 for all the elements, and make p = 2 for the 7th element, the solution is still correct. However, if I use p=2 for all the elements and make p=3 for the 7th element, the solution suddenly becomes completely wrong. Choosing the 7th element is just arbitrary, other choices also show the same problem.

Since it works for the combination of p=1 and 2, I do not think hanging nodes are a problem. I looked at the hp-paper as well and do not see something that might be going wrong related to the hanging nodes. I am still trying to narrow down the problem and come up with a simplified example, however, I still thought of raising this question. Based on the above scenario, can someone figure out what might be going wrong?

Best regards
Deepak

[Wolfgang Bangerth]

On 10/05/2016 09:03 AM, Deepak Gupta wrote:
>
>
> Since it works for the combination of p=1 and 2, I do not think hanging nodes
> are a problem. I looked at the hp-paper as well and do not see something that
> might be going wrong related to the hanging nodes. I am still trying to narrow
> down the problem and come up with a simplified example, however, I still
> thought of raising this question. Based on the above scenario, can someone
> figure out what might be going wrong?

I am now aware of anything that would cause this, but it certainly sounds like
a bug. A small testcase would definitely be useful!

In the past, I have debugged problems such as this by interpolating a known
and simple function onto the mesh, and then applying the constraints via
ConstraintMatrix::distribute() to the vector with the interpolant. If the
constraints are wrong, you will see this if you just output this as a field
(with sufficiently high n_subdivisions passed to DataOut::build_patches()).

Best
W.

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

[Deepak Gupta]

Thanks Wolfgang for the reply. I created a small scalar example by modifying Step-6 and to my surprise, it works very well (can send if someone needs it). So, I have not figured out the error yet. I am still trying to narrow down the problem.

By the way, I am constructing the cell matrices using B'DB, where B is the strain-displacement matrix and D is the constitutive relation. Then I add it to the system matrix using distribute_local_to_gobal function. Since B and D are the same for global and adaptive case, I do not expect any difference in the cell matrices. The system_rhs is correct since I have only a point load and I can check it by printing the vector. Thus, the possible issue can be only in the assembly (for the local refinement one). My problem is vector-values. In case this information can lead to any other suggestion from someone, it would be great.

Meanwhile I am going to extend the simple example for a basic version of the elastic problem I am trying to solve. Hope then I can figure out the error.

Best regards

Deepak

--
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+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Wolfgang Bangerth]

Deepak,

> By the way, I am constructing the cell matrices using B'DB, where B is
> the strain-displacement matrix and D is the constitutive relation. Then
> I add it to the system matrix using distribute_local_to_gobal function.
> Since B and D are the same for global and adaptive case, I do not expect
> any difference in the cell matrices. The system_rhs is correct since I
> have only a point load and I can check it by printing the vector. Thus,
> the possible issue can be only in the assembly (for the local refinement
> one). My problem is vector-values. In case this information can lead to
> any other suggestion from someone, it would be great.

I don't have any experience with this approach, but at the end of the
day, you get a local matrix A_K = B_K' D_K B_K for each cell K, and
that should be all that matters. I assume you mean
ConstraintMatrix::distribute_local_to_gobal, not
cell->distribute_local_to_gobal, right?

> Meanwhile I am going to extend the simple example for a basic version of
> the elastic problem I am trying to solve. Hope then I can figure out the
> error.

Good plan.

[Deepak Gupta]

Hi Wolfgang,

Yes it is ConstraintMatrix::distribute_local_to_global.

I have finally come up with a small example that clearly shows the issue I am facing. Attached is the code and 3 outputs. I have a rectangle domain of 4X8 elements. A point load (0, 1) is applied at the second vertex from the top at the right boundary and the left boundary is fixed. The y-displacement will be maximum at the load point. Now I look at 3 cases:

1) Using globally biquadratic elements -> displacement value is 420.2 (as in biquadratic.png).
2) Using globally bicubic elements -> displacement value is  426.8 (as in bicubic.png).

3) Using biquadratic for all elements, except that the element with index 6 is bicubic, displacement value is 61.45 (as in mix.png).

I would expect the result of 3) to be between 1) and 2) ideally. I have tested it for finer meshes, but such issues exists. Could someone help me in locating the problem. Thanks in advance.

Best regards

Deepak

[Jean-Paul Pelteret]

Hi Deepak,

You've chosen a difficult example for us to evaluate, because your implementation of the plane stress condition is rather hard to follow. Can you provide a reference from which you're implementing this?

At the very least your implementation of "add_point_source_to_rhs" looks suspicious to me. The ordering of the vertices has, in general, nothing to do with the ordering of the DoFs. You should use the vertex_dof_index function to retrieve the DoF associated with a vertex. There are a number of posts that discuss the use of this function.

J-P

[Deepak Gupta]

Dear JP,

Thanks for looking at the example. The code itself is not a simple example that I implemented. Rather it is a part of my exact problem and I have cut it out from the rest of the code to demonstrate my problem. I am working on design optimization and this part goes in there as the modeling part.

Indeed the compilation of the B-matrix might look complicated. However, a closer look might reveal that B matrix is compiled as shown in http://homepages.rpi.edu/~des/4NodeQuad.pdf, page 2, slide 3. This link is just a  random example I took. I hope there is no error in  that part of the code since I have been using it and it worked in the past. Also, if that formulation would be wrong, I would wonder why the globally fine mesh would give the correct results. The values 420, 428 etc are correct since they are very close to a MATLAB version of FEM code that I compared with.

Thanks for pointing out the vertex_dof_index function. However, I have some question related to it now:

1) When I describe the point load using my approach, I get the load where I wanted (see my_approach_point_load.png).

2) However, when I use vertex_dof_index function, the load gets attached to a different component of a wrong vertex (see new_point_load.png). I am wondering how I should use this function, or may be there is something else wrong with my code that makes it wrong.

I used following lines in the code:
                system_rhs(cell->vertex_dof_index(v, 0)) = 0;
                system_rhs(cell->vertex_dof_index(v, 1)) = 1;

3) In addition, I would like to know whether I can use

Ve

ctorTools::create_point_source_vector for the same? This function used to help me for bilinear elements.

Best regards

Deepak

[Jean-Paul Pelteret]

Dear Deepak,

Its always easier to help debug a well document, clear and concise example. It was less about checking the implementation than understanding what the code is doing and how all of the parts interact. After all, we can only work with what you give us.

Anyway, I've noticed that you call "add_point_source_to_rhs" after you've applied the hanging node constraints and boundary conditions via hanging_node_constraints.distribute_local_to_global. This might be the problem because the matrix and RHS contribution is modified when the local constraints are applied. I would suggest that you implement this as a standard traction boundary condition and then add its local contribution to the global RHS vector as per usual. This would eliminate another possible source of error.

Yes, as best I can discern VectorTools::create_point_source_vector should do the same job for you (i.e. create the equivalent RHS contribution as a point load). I've never used it though so I don't know what the implications of using it in this context are (i.e. the interaction of hp-constraints).

Regards,

J-P

On Friday, October 7, 2016 at 8:41:42 AM UTC+2, Deepak Gupta wrote:

Dear JP,

Thanks for looking at the example. The code itself is not a simple example that I implemented. Rather it is a part of my exact problem and I have cut it out from the rest of the code to demonstrate my problem. I am working on design optimization and this part goes in there as the modeling part.

Indeed the compilation of the B-matrix might look complicated. However, a closer look might reveal that B matrix is compiled as shown in http://homepages.rpi.edu/~des/4NodeQuad.pdf, page 2, slide 3. This link is just a  random example I took. I hope there is no error in  that part of the code since I have been using it and it worked in the past. Also, if that formulation would be wrong, I would wonder why the globally fine mesh would give the correct results. The values 420, 428 etc are correct since they are very close to a MATLAB version of FEM code that I compared with.

Thanks for pointing out the vertex_dof_index function. However, I have some question related to it now:

1) When I describe the point load using my approach, I get the load where I wanted (see my_approach_point_load.png).

2) However, when I use vertex_dof_index function, the load gets attached to a different component of a wrong vertex (see new_point_load.png). I am wondering how I should use this function, or may be there is something else wrong with my code that makes it wrong.

I used following lines in the code:
                system_rhs(cell->vertex_dof_index(v, 0)) = 0;
                system_rhs(cell->vertex_dof_index(v, 1)) = 1;

3) In addition, I would like to know whether I can use

Ve

ctorTools::create_point_source_vector for the same? This function used to help me for bilinear elements.

Best regards

Deepak

[Deepak Gupta]

Thanks JP. Indeed a concise example should help and that is where I started. I took a scalar-valued example, however, it did not show the same problem as this one. Thus, I decided to move towards this case. I take this blame that I should have an even simpler case, and I hope I can come up soon.

In the meantime, I will take the suggestions into account and see if the issue goes away.

"Anyway, I've noticed that you call "add_point_source_to_rhs" after you've applied the hanging node constraints and boundary conditions via hanging_node_constraints.distribute_local_to_global. This might be the problem because the matrix and RHS contribution is modified when the local constraints are applied. I would suggest that you implement this as a standard traction boundary condition and then add its local contribution to the global RHS vector as per usual. This would eliminate another possible source of error."

So, would it be correct, if I move my add_point_source_to_rhs function above create_local_to_global? Also, to compare with the benchmarks, I need to keep it a point load. Would that still be possible with traction boundary condition?

Best,

Deepak

[Jean-Paul Pelteret]

Ok, so its a benchmark test not one you've come up with yourself. Fair enough...

What I'm suggesting is that you implement the RHS contribution from the load within the local assembly loop. You can have a look at how the point source is implemented in VectorTools::create_point_source_vector and duplicate this within your assembly routine (you could probably either do this localised to one cell or distribute the load contribution over the neighbouring cells). This way when you call distribute_local_to_global this contribution gets scaled correctly. Why you can't do this with your current function "add_point_source_to_rhs" is because it directly modifies the global system RHS vector.

I hope that this explains better what I think you might need to do to correct the issue.

J-P

[Deepak Gupta]

Thanks JP for the elaborate clarification. I quickly tried to test a distributed source term commenting out the point load function and adding the following before distribute_local_to_global:
      
 //Calculating cell_rhs
        for(unsigned int i = 0; i < dofs_per_cell; ++i){
            const unsigned int component_i = cell->get_fe().system_to_component_index(i).first;
            for(unsigned int q_point = 0; q_point < n_q_points; ++q_point){
                double body_load = 0;
                if (component_i == 1 && cell_itr%4 == 3) body_load = 0.5;
                cell_rhs(i) += fe_values.shape_value(i, q_point) *
                        body_load * //rhs_values[q_point](component_i) *
                        fe_values.JxW(q_point);
            }
        }

However, the results are still similar, although magnitudes differ (globally_p_2.png and locally_p_3.png). The cases with global and local refinement are attached. The strange thing is that the symmetry in absolute value of x-displacement is lost when one element is of higher order. I would expect a bit of difference but not this much.

Do you think the source term implementation is correct? I added force term only in the two cells on the ride side.

Best regards

Deepak

[Jean-Paul Pelteret]

Hi Deepak,

Its not pretty but at a glance this does look plausibly correct. There are a few more things that you should check:

- Add the boundary conditions directly to the constraints matrix, instead of using MatrixTools::apply_boundary_values. You can just use a ZeroFunction on boundary 42 instead of the function you've defined here.

- You've created a temporary quadrature object "QGauss<2> temp_quad(cell->active_fe_index()+3);" which is suspicious. Remove this. You can retrieve the active quadrature formula from FEValues.

- Run this in debug mode to check that there are no other mistakes that we can't pick up by just looking at the code.

If you do all of this and still have no success, then post the updated code and I'll fiddle with it later.

J-P

[Deepak Gupta]

Thanks JP,

I took the above points into my code except that I do not know/understand how to add boundary conditions in the constraint matrix without using the apply_boundary_values. Could you divert me to an example or discussion where this is done.

Best

Deepak

[Jean-Paul Pelteret]

See step-27:

VectorTools::interpolate_boundary_values (dof_handler,
                                            0,
                                            ZeroFunction<dim>(),
                                            constraints);

[Deepak Gupta]

Thanks.

Just some add on, I figured out that the replacement for the lines commented below are the ones in bold. Here it does not make sense to me why v-1, 1 and 2 work, rather than v, 0 and 1 (as seen in google-group discussion). May be this information is of use to you in back-tracing some error in my code.

By replacement, I mean same system_rhs vector, as viewed in paraview.

        for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v){
            if (fabs(cell->vertex(v)[0] - 2.0) < 1e-12 && fabs(cell->vertex(v)[1] - 0.5) < 1e-12){
                system_rhs(cell->vertex_dof_index(v-1, 1)) = 0;
                system_rhs(cell->vertex_dof_index(v-1, 2)) = 1;
                //system_rhs[global_dof_indices[2*v]] = 0;
                //system_rhs[global_dof_indices[2*v + 1]] = 1;
            }
        }

I will check and get back with the updated code, in case I do not figure out the error :)

Best regards

Deepak

[Deepak Gupta]

Dear JP and Wolfgang,

Thanks for your help. I could finally figure out the bug. Credit to JP who doubted by B matrix formulation and indeed it went wrong due to the DOF numbering system in dealii.

To be more precise, I had assumed that all dofs associated with a certain vertex are arranged together (one after the other) in the  whole dof_list. However, this was not true beyond p=2 and that is where my code failed.

Best regards

Deepak

[Jean-Paul Pelteret]

Great! I'm glad that you found out the problem! 

Unfortunately this is the sort of issue that trips up many "novice" users of deal.II. We recognise that its hard to get into the habit of writing code that is completely free of assumptions of this sort, since its often quicker (but not necessarily easier) to hardcode in the patterns that one notices for any given dimension, FE type etc. However, we strongly believe that its worth investing the few extra moments to structure the code correctly up-front because you save time in the long run when your code gets more complex. You'll notice that in the documentation and the lectures, the concept of authoring general code is very strongly reinforced. This is especially relevant as many authors contribute to many different parts of the code, and we each have our own philosophy as to what a reasonable ordering is for whatever task is thats being implemented (or one approach is more practically implementable than another). Take for instance FE_Q vs FE_DGQ: they're so similar in nature but very different in implementation. Much effort has been taken to provide helper and compatibility functions to work around this issue. By using these functions wherein you can query local DoF components, FE groups etc., any troubles related to this is completely avoided.

I recommend that anyone learning deal.II take heed of this example, and consider investing a little more time in the best quality (read "most defensive") code possible rather than one just good enough to achieve the current task.

J-P

[Jean-Paul Pelteret]

Oh, and I should say "kudos" to you for sticking with it and finding the bug yourself. I like to think that its all part of the learning experience :-)