idaholab / moose

Multiphysics Object Oriented Simulation Environment
https://www.mooseframework.org
GNU Lesser General Public License v2.1
1.79k stars 1.05k forks source link

Computing material properties on displaced mesh causes segfault #13151

Open reverendbedford opened 5 years ago

reverendbedford commented 5 years ago

Bug Description

Computing material properties on the displaced mesh leads to a segfault when the material is evaluated on a boundary. Top of the stack trace is:

#0  0x00007ffff6cf8240 in std::__cxx1998::vector<libMesh::Point, std::allocator<libMesh::Point> >::end() const (this=0x20)
    at /usr/include/c++/4.8.2/bits/stl_vector.h:566
#1  0x00007ffff6da7384 in std::__cxx1998::vector<libMesh::Point, std::allocator<libMesh::Point> >::empty() const (this=0x20)
    at /usr/include/c++/4.8.2/bits/stl_vector.h:735
#2  0x00007ffff6da6d44 in libMesh::QBase::n_points() const (this=0x0)
    at /home/messner/projects/moose/scripts/../libmesh/installed/include/libmesh/quadrature.h:129
#3  0x00007ffff409f6fe in Material::computeProperties() (this=0xa08128)
...

which seems to imply that it's a problem accessing the list of quadrature points, maybe in libmesh?

Steps to Reproduce

  1. Compile tensor_mechanics
  2. Run the following input
# Simple plane strain test

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Variables]
      [./disp_x]
      [../]
      [./disp_y]
      [../]
[]

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 4
  ny = 4
[]

[Modules/TensorMechanics/Master]
  [./all]
    strain = FINITE
    planar_formulation = PLANE_STRAIN
    add_variables = true
    generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy strain_zz'
  [../]
[]

[Functions]
  [./pull]
    type = ParsedFunction
    value ='100 * t'
  [../]
[]

[BCs]
  [./leftx]
    type = PresetBC
    boundary = right
    variable = disp_x
    value = 0.0
  [../]
  [./bottomy]
    type = PresetBC
    boundary = bottom
    variable = disp_y
    value = 0.0
  [../]
  [./pull]
    type = FunctionNeumannBC
    boundary = top
    variable = disp_y
    function = pull
  [../]
[]

[Materials]
  [./elastic_stress]
    type = ComputeFiniteStrainElasticStress
  [../]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = 0.33
    youngs_modulus = 1000.0
  [../]
  [./dummy]
      type = GenericConstantMaterial
      prop_names = 'dummy'
      prop_values = '1'
      use_displaced_mesh = true
  []
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient

  solve_type = 'newton'
  line_search = none

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  l_max_its = 3
  l_tol = 1e-15
  nl_max_its = 30
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-8

  start_time = 0.0
  dt = 1.0
  dtmin = 1.0
  end_time = 5.0
[]

[Outputs]
  exodus = true
[]

To work around either:

  1. Remove use_displaced_mesh = true from the dummy material
  2. Add if (_bnd) return; to Material::computeProperties(), which is how I know it has to do with evaluating on boundaries.

Impact

I can use the if (_bnd) return workaround in my app, but I'm a little worried about what that implies for evaluating boundary conditions.

permcody commented 5 years ago

Thanks for the detailed report and test case!