Step 44: imported 2d mesh fails to converge to solution

[James Gorman]

To Whom It May Concern,

I am learning deal.ii to work on problems in solid mechanics. In the interest of computational efficiency, I plan to run most problems in 2D (plane stress assumption). Since the interesting problems I am concerned about require a non-trivial mesh, I learned how to use Gmsh based on the suggestion of Step-49.

With that in mind, I ran Step-44 in 2D and 3D (no change to the tutorials) and got results. I then used gmsh to generate a hypercube mesh in 2D and 3D, modifying Step-44 using the Step-5 as a guide. While the 3D calculation with the imported mesh successfully converged, the 2D calculation with imported mesh did not converge. The 2D calculation failed with the residual from 10^5-10^16 (I tried multiple initial mesh densities).

Looking at the error output, it states that when a residual is large the matrix may be invertible or the wrong solver may be used. Since the only thing I changed was the origin of the mesh, I don't think either of these explanations should be the case.

As a sanity check, I ran a ~5 of the other tutorials with imported meshes. None of the other tutorials I tried ran into a similar issue when I imported a mesh from Gmsh. The issue also stuck around whether or not I used a structured or unstructured mesh, or changed the shape from a square to a rectangle.

Has anyone else ran into this issue? If so, how did you fix it?

Sincerely,

James Gorman

[Jean-Paul Pelteret]

Dear James,

To me, the most likely point of failure here relates to the application of boundary conditions. It might be that one of the constrained boundaries is not traction free, and so the problem is not free of rigid body modes. Can you please confirm that you apply the same numbering scheme to the boundary in 2d as is used in the tutorial? You can confirm the setup in either case by looping over all boundary faces after the mesh is created, and printing their ID. I think that this is the best starting point to diagnose what the problem could be. If that’s not the problem, then can you please provide an example mesh that fails, and confirm which version of deal.II you’re using.

Best,

Jean-Paul

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/f8bd877d-5d1c-477d-8d84-70903f15efb2%40googlegroups.com.

[Jean-Paul Pelteret]

It might be that one of the constrained boundaries is not traction free...

Correction: This should have read "It might be that one of the constrained boundaries is now traction free..."

[James Gorman]

Jean-Paul,

The boundary ID output from the four different calculations are as follows (before any boundary IDs are set to 6):

1) 3D, no imported mesh: equal numbers for boundary IDs 0, 1, 2, 3, 4, and 5. This calculation converges to a solution (as expected). 

2) 2D, no imported mesh: equal numbers for boundary IDs 0, 1, 2, and 3. This calculation converges to a solution (as expected). 

3) 3D, imported mesh: some boundary IDs are 1 and 2, while most boundary IDs are 3. This calculation converges to a solution.

4) 2D, imported mesh: all boundary IDs are 0. This calculation does not converge to a solution.

After cycling through all of the face cells and using set_boundary_id(), each calculation (correctly) has 6 as a boundary ID for some number of faces.

With the information above and your comments, I take it that boundary ID = 0 means the boundary is traction-free and no dirichlet BC is specified? Assuming this is the case, then it makes sense there would be no convergence. It also means that the 3D calculation with the imported mesh is not the same as the original. I have not yet determined why I cannot read the .vtu output files in Paraview, which is why I did not know this before (I could only tell it converged to a solution).

I do find it interesting that the two meshes were not imported with the same boundary IDs. I guess in the future I should assume that all boundary IDs are zero and reset them for each calculation? Also, is there a quick reference guide for what each boundary ID means? I have not yet found such a reference on the deal.ii website.

Thank you for your help,

James

[Daniel Arndt]

James,

I do find it interesting that the two meshes were not imported with the same boundary IDs. I guess in the future I should assume that all boundary IDs are zero and reset them for each calculation? Also, is there a quick reference guide for what each boundary ID means? I have not yet found such a reference on the deal.ii website.

There is no global meaning for any boundary id. Every code can decide what to do for faces with a specific boundary id. Have a look, e.g., where interpolate_boundary_values is called and have a look at the first tutorial steps, in particular step-3.

Best,

Daniel

[James Gorman]

Daniel,

Thank you for the clarification. The modified tutorial is working now.

James