Source of a constant residual error when a triangulation has hanging nodes

[amir.pa...@gmail.com]

Dear all

I have a code that works fine. Now, I am working to change the code structure to use “workstream”. The new code works fine without adaptive refinement. However, when the mesh is refined and the triangulation has hanging node, a constant residual error remains during Newton iteration. I cant find source of the problem.

There are two “constraints matrix”. On of them for “Newton update” and one for “hanging nodes”. The constraints matrix are merged and is used to determine error as:

BlockVector<double> error_ud(dofs_per_block);

for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)

if (!constraints_update.is_constrained(i))

error_ud(i) = system_rhs(i);

error_update.norm = error_ud.l2_norm();

Regards

Pasha

[Daniel Arndt]

Pasha,

I have a code that works fine. Now, I am working to change the code structure to use “workstream”. The new code works fine without adaptive refinement. However, when the mesh is refined and the triangulation has hanging node, a constant residual error remains during Newton iteration. I cant find source of the problem

You are providing way too few information to expect useful answers. What are you doing in your code?
For which part in your code do you want to change to WorkStream?
What exactly did you change? How is this related to whether you are doing adaptive refinement or not?

Best,
Daniel

[amir.pa...@gmail.com]

Dear Daniel

The code solve a solid mechanic model with a plastic material model. The WorkStream is used for matrix and residual assembly (like step-44). So I just add some new classes  (PerTaskData and ScratchData) and the most parts of the code is not changed. The code works similar to the its serial version without any problem when the local adaptive refinement is not active. During a load increment, let the solution to time step Tn be given. The full problem is solved. First, I refine the mesh based on the new solution and interpolate the old solution (at Tn ) onto the new mesh. Then I solve for the solution at Tn+1 again, but on the refined mesh. This process is tested in the previous version of the code. However, in the new code after the first refinement, the solution for Tn+ has a constant residual error during Newton iteration. 

Thank you

P.

[amir.pa...@gmail.com]

Dear Denis

I checked the the system residual assembly process. For a triangulation with 28 cells, the dof index 17 is belong to cells 2, 5, 8 and 27. The corresponding cell_rhs values for this index are,

cell        cell_rhs value

2            2.88462e+00

5            2.88462e+00

8           -2.88462e+00

27         -1.44231e+00

This dof is not constrained. When the following assembly method is used:

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

system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i);

the system_rhs(17) = 1.4423076, however, using the following method,

constraints_update.distribute_local_to_global(data.cell_rhs,

data.local_dof_indices, system_rhs);

the system_rhs(17) = -1.33e-15. This is the source of the problem because I have used the second method in the previous code and the first one (similar to step-44) in the new code. I cant understand why the results are different.

Regards

P.

On Wednesday, February 1, 2017 at 8:24:51 PM UTC+3:30, Daniel Arndt wrote:

[Daniel Arndt]

Pasha,

 

I checked the the system residual assembly process. For a triangulation with 28 cells, the dof index 17 is belong to cells 2, 5, 8 and 27. The corresponding cell_rhs values for this index are,

cell        cell_rhs value

2            2.88462e+00

5            2.88462e+00

8           -2.88462e+00

27         -1.44231e+00

This dof is not constrained. When the following assembly method is used:

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

system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i);

the system_rhs(17) = 1.4423076, however, using the following method,

constraints_update.distribute_local_to_global(data.cell_rhs,

data.local_dof_indices, system_rhs);

the system_rhs(17) = -1.33e-15. This is the source of the problem because I have used the second method in the previous code and the first one (similar to step-44) in the new code. I cant understand why the results are different.

This looks like DoF 17 is constrained, possibly by hanging node constraints. Can you recheck the constraints stored in `constraints_update` via `constraints_update.print(std::cout)`?
In any case, you need to use a ConstraintMatrix for assembling the right-hand side either via calling `distribute_local_to_global` on each cell or calling `condense` after the loop over all cells. After solving you would call `constraints.distribute`.
Otherwise, you would not take the hanging node constraints into account.
Are you missing the call to `constraints.condense(system_matrix, system_rhs);` in the first approach?

Best,
Daniel

[amir.pa...@gmail.com]

Dear Daniel

After solving constraints_update.distribute is called. In addition, after assembling the system matrix and residual vector using the first method, constraints_update.condense is called before solving the system. After that the solution is updated and the new residual is calculated based on the first method, so in error calculation the constraints are considered as,

BlockVector<double> error_res(dofs_per_block);

for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)

if (!constraints_update.is_constrained(i))

error_res(i) = system_rhs(i);

error_residual.norm = error_res.l2_norm();

constraints update

1 = 0

3 = 0

11 = 0

21 = 0

22 10: 5.00000e-01

22 14: 5.00000e-01

23 15: 5.00000e-01

24 14: 5.00000e-01

24 16: 5.00000e-01

25 15: 5.00000e-01

25 17: 5.00000e-01

26 16: 5.00000e-01

26 40: 5.00000e-01

27 17: 5.00000e-01

27 41: 5.00000e-01

42 = 0

43 = 1.00000e-04

45 = 1.00000e-04

46 = 0

49 = 1.00000e-04

51 = 1.00000e-04

53 = 1.00000e-04

55 = 0

56 40: 5.00000e-01

57 41: 5.00000e-01

57 59: 5.00000e-01

58 = 0

69 = 0

70 = 0

72 = 0

78 = 0

79 = 0

80 = 0

93 87: 5.00000e-01

93 89: 5.00000e-01

94 89: 5.00000e-01

94 90: 5.00000e-01

95 90: 5.00000e-01

95 102: 5.00000e-01

110 102: 5.00000e-01

110 111: 5.00000e-01

constraints hanging nodes

22 10: 5.00000e-01

22 14: 5.00000e-01

23 11: 5.00000e-01

23 15: 5.00000e-01

24 14: 5.00000e-01

24 16: 5.00000e-01

25 15: 5.00000e-01

25 17: 5.00000e-01

26 16: 5.00000e-01

26 40: 5.00000e-01

27 17: 5.00000e-01

27 41: 5.00000e-01

56 40: 5.00000e-01

56 58: 5.00000e-01

57 41: 5.00000e-01

57 59: 5.00000e-01

93 87: 5.00000e-01

93 89: 5.00000e-01

94 89: 5.00000e-01

94 90: 5.00000e-01

95 90: 5.00000e-01

95 102: 5.00000e-01

110 102: 5.00000e-01

110 111: 5.00000e-01

Thank you.

Pasha