FileMeshGenerator and uniform_refine produce weird results

[maxnezdyur]

Bug Description

I had some weird results with an exodus file of mine when refining it. As I increased uniform_refine I got weird results. I found MWE using an exodus file MOOSE already had. Seems like child elements are not C0?

Steps to Reproduce

Exodus file exists in MOOSE at modules/solid_mechanics/test/tests/volumetric_locking_verification/42_node_side.e

# 2D cook's membrane problem with a trapezoid
# that is fixed at one end and is sheared at
# other end. Poisson's ratio is 0.4999.

# Using Quad4 elements and no volumetric locking,
# vertical displacement at top right corner is 3.78
# due to locking.

# Using Quad4 elements with volumetric locking, vertical
# dispalcement at top right corner is 7.78.

# Results match with Nakshatrala et al., Comp. Mech., 41, 2008.

# Check volumetric locking correction documentation for
# more details about this problem.

[GlobalParams]
  displacements = 'disp_x disp_y'
  # volumetric_locking_correction = true
[]

[Mesh]
  [generate]
    type = FileMeshGenerator
    file = 42_node_side.e
  []
  second_order = true
  uniform_refine = 3
[]

[Variables]
  [disp_x]
    order = SECOND
  []
  [disp_y]
    order = SECOND
  []
[]
[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        incremental = true
      []
    []
  []
[]

[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = 1
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  [Nuem]
    type = NeumannBC
    variable = disp_y
    boundary = 2
    value = 6.25
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 250.0
    poissons_ratio = 0.4999
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = 'NEWTON'

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  num_steps = 1
  nl_abs_tol = 1e-5
[]

[Postprocessors]
  [a_disp_y]
    type = PointValue
    variable = disp_y
    point = '48.0 60.0 0.0'
  []
[]

[Outputs]
  exodus = true
  file_base = test_higher
[]

Impact

Mesh refinement not correct.

[Optional] Diagnostics

No response

[lindsayad]

are you confident that this is not a visualization issue?

[maxnezdyur]

Pretty confident. The results of this problem should converge to something like 0.8 for the post processor. The results are like 0.7, 0.75, then 1.8, as I increase the refinement.

[lindsayad]

There seems to be a bad interaction between uniform_refine and second_order = true. Running in devel mode I get the backtrace

(gdb) bt
#0  0x00007ffff1a4ee90 in libMesh::MacroFunctions::report_error(char const*, int, char const*, char const*, std::basic_ostream<char, std::char_traits<char> >&)@plt ()
   from /home/lindad/projects/moose-no-cuda/scripts/../libmesh/installed/lib/libmesh_devel.so.0
#1  0x00007ffff2605c75 in libMesh::MeshRefinement::add_node (this=0x7fffffffc510, parent=..., child=<optimized out>, 
    node=5, proc_id=0) at ../src/mesh/mesh_refinement.C:157
#2  0x00007ffff229c993 in libMesh::Elem::refine (this=0x5555559c6d00, mesh_refinement=...)
    at /usr/include/c++/12/bits/unique_ptr.h:191
#3  0x00007ffff260ab02 in libMesh::MeshRefinement::_refine_elements (this=this@entry=0x7fffffffc510)
    at ../src/mesh/mesh_refinement.C:1556
#4  0x00007ffff260b30e in libMesh::MeshRefinement::uniformly_refine (this=this@entry=0x7fffffffc510, n=n@entry=3)
    at ../src/mesh/mesh_refinement.C:1697
#5  0x00007ffff5900550 in Adaptivity::uniformRefine (mesh=0x5555559ce290, level=3, level@entry=4294967295)
    at /home/lindad/projects/moose-no-cuda/framework/src/base/Adaptivity.C:289
#6  0x00007ffff634da5d in SetupMeshCompleteAction::act (this=0x5555557da2c0)
    at /home/lindad/projects/moose-no-cuda/framework/src/actions/SetupMeshCompleteAction.C:78
#7  0x00007ffff634b460 in Action::timedAct (this=0x5555557da2c0)
    at /home/lindad/projects/moose-no-cuda/framework/src/actions/Action.C:87
#8  0x00007ffff6359b9d in ActionWarehouse::executeActionsWithAction (this=this@entry=0x5555557feee8, 
    task="uniform_refine_mesh") at /home/lindad/projects/moose-no-cuda/framework/src/actions/ActionWarehouse.C:387
#9  0x00007ffff6359ddf in ActionWarehouse::executeAllActions (this=this@entry=0x5555557feee8)
    at /home/lindad/projects/moose-no-cuda/framework/src/actions/ActionWarehouse.C:346
#10 0x00007ffff58851a8 in MooseApp::runInputFile (this=0x5555557fe7d0)
    at /home/lindad/projects/moose-no-cuda/framework/src/base/MooseApp.C:1105
#11 0x00007ffff5889caa in MooseApp::run (this=0x5555557fe7d0)
    at /home/lindad/projects/moose-no-cuda/framework/src/base/MooseApp.C:1534
#12 0x00005555555573f6 in Moose::main<MooseTestApp> (argc=4, argv=0x7fffffffe708)
    at /home/lindad/projects/moose-no-cuda/framework/build/header_symlinks/MooseMain.h:47
#13 0x00005555555574bd in main (argc=<optimized out>, argv=<optimized out>)
    at /home/lindad/projects/moose-no-cuda/test/src/main.C:16

@roystgnr is this something you'd want to look into? If not, I'll investigate

[roystgnr]

This is really astonishing. "Refine a second-order mesh" is a really common thing to do. Let me see if I can replicate.

[lindsayad]

I wouldn't be surprised if it's something dumb we're doing in MOOSE

[roystgnr]

Even in MOOSE it seemed like "go second order, refine" would be both straightforward (you call one libMesh function, then another, and order doesn't matter) and frequently used. I figured the problem had to be something wrong with that particular Exodus file, but Paraview makes it look both perfectly well-formed and too simple for anything to be subtly wrong.

I can replicate it and I've got it in the debugger now, and the first thing that screams at me is that I'm seeing the error message from a QUADSHELL8!? I was expecting to see a simple QUAD9. I don't think we have much test coverage of the "shell" elements so there could totally be an issue there.

[roystgnr]

We do build a QUADSHELL4 if the ExodusII file specifies a SHELL4 rather than a QUAD4, and that would p-refine to a QUADSHELL8.

[roystgnr]

Yeah, ncdump 42_node_side.e | grep type sees connect1:elem_type = "SHELL4" ;

[roystgnr]

I'd say the workaround here is to use Quad8 or Quad9 instead, if that's appropriate? But I use shell elements so little that we don't even have a QUADSHELL4->QUAD4 conversion utility, and based on that file's header (which seems to have shellface sidesets?) I'm guessing the answer to "if that's appropriate" is "it isn't"?

The fix is going to take me a little while, but I'm guessing this is a libMesh bug. Man, I don't even have INSTANTIATE_ELEMTEST(QUADSHELL4) or INSTANTIATE_ELEMTEST(QUADSHELL8) in our unit tests; it's like I'd completely forgotten those existed.

[maxnezdyur]

I think cubit defaults to QUADSHELL4 if you don't specify an element type, which is why I ran into the problem.

[lindsayad]

In that case, it's kind of shocking we haven't run into this before

[roystgnr]

Adding those missing unit tests, it looks like refinement on QUADSHELL4 is working; only QUADSHELL8 is failing. So it really is only the combination of shell inputs and all_second_order and` refinement that breaks.

[roystgnr]

What I can't figure out yet is why this doesn't affect QUAD8. The point of failure is somewhere that I'd expect to hit QUAD8 too, and Quad8Shell is just a Quad8 subclass that changes the returned type().

[roystgnr]

Okay, I get it. Sylvain added Quad8Shell but no Quad9Shell, so whereas Quad8 "cheats" at refinement by looking at the Quad9 center node to establish bracketing_nodes connections, Quad8Shell just fails there.

[roystgnr]

I'm just going to add Quad9Shell. It's got an ExodusII type, SHELL9, so maybe people will find it directly useful too.

[roystgnr]

This file seems to be refining fine with https://github.com/libMesh/libmesh/pull/3886 for me. Hopefully we'll get that into master in time for a MOOSE libMesh submodule update this week.

[roystgnr]

Gah, that "Hopefully this will fix" wasn't intended as a keyword. This isn't fixed until we get the libMesh submodule updated.

[roystgnr]

This should be fixed as of #28014