How to ensure the displaced mesh is used

[keassyguang]

Hello MOOSE Community,

I'm developing a custom material class that inherits from ADMaterial. The coordinates of the current quadrature point will be used as a condition to assign a specific value to some created variables. What I want to know is whether the coordinates of _q_point[_qp](0), _q_point[_qp](1) and _q_point[_qp](2) return the displaced mesh or the undisplaced mesh. Since I need the coordinates of the quadrature point under the displaced mesh, what should I do to ensure that the displaced mesh is always used?

[GiudGiud]

For many objects, all you need to do is setting use_displaced_mesh = true

Is that an option for your case?

[keassyguang]

These inputs dont run without using your app. I ll have to use a machine on which I can install an external app. I ll be delayed

Thanks a lot.

I have created another non-working example that doesn't require my app. This example can only run the first two or three steps, but I think it is enough to illustrate the problem as the adaptivity does not work and the calculated gradient jump is also close to zero with the displaced mesh.

Key results:

1. with use_displaced_mesh = true, the adaptivity does not work, gradient jump value close to zero;
2. comment out use_displaced_mesh = true, the adaptivity works but is aborted with a converge issue.

Cannot take square root of a number less than or equal to zero in the calculation of C3_test for the Rashid approximation for the rotation tensor. This zero or negative number may occur when elements become heavily distorted.
To recover, the solution will fail and then be re-attempted with a reduced time step.

  Nonlinear solve did not converge due to DIVERGED_FNORM_NAN iterations 3
 Solve Did NOT Converge!
Aborting as solve did not converge

Input file

[Mesh]
  type = GeneratedMesh
  displacements = 'disp_x disp_y'
  dim = 2
  xmin = 0
  ymin = 0
  xmax = 2
  ymax = 1
  nx = 10
  ny = 5
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
  use_displaced_mesh = true
[]

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

[AuxVariables]
  [von_mises]
    order = FIRST
    family = MONOMIAL
  []
[]

[AuxKernels]
  [von_mises]
    type = ADRankTwoScalarAux
    variable = von_mises
    rank_two_tensor = stress
    scalar_type = VonMisesStress
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Kernels]
  [./stress_x]
    type = ADStressDivergenceTensors
    displacements = 'disp_x disp_y'
    component = 0
    variable = disp_x
  [../]
  [./stress_y]
    type = ADStressDivergenceTensors
    displacements = 'disp_x disp_y'
    component = 1
    variable = disp_y
  [../]
[]

[BCs]
  [./fix_bottom_x]
    type = ADDirichletBC
    variable = disp_x
    boundary = bottom
    value = 0
  [../]
  [./fix_bottom_y]
    type = ADDirichletBC
    variable = disp_y
    boundary = bottom
    value = 0
  [../]
  [top_pull_y]
    type = FunctionDirichletBC
    variable = disp_y
    function = 't'
    boundary = 'top'
  []
[]

[Materials]
  [./elasticity]
    type = ADComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 10e9
  [../]
  [./strain]
    type = ADComputeFiniteStrain
    displacements = 'disp_x disp_y'
  [../]
  [./stress]
    type = ADComputeFiniteStrainElasticStress
  [../]
[]
  
[Executioner]
  type = Transient
  dt = 0.05
  solve_type = 'PJFNK'

  petsc_options_iname = -pc_hypre_type
  petsc_options_value = boomeramg

  dtmin = 0.01
  num_steps = 10
[]

[Outputs]
  exodus = true
[]

[Adaptivity]
  max_h_level = 4
  initial_marker = combo
  marker = combo
  
  [Indicators]
    [jump_von_mises]
      type = GradientJumpIndicator
      variable = von_mises
    []
  []

  [Markers]
    [errorfrac]
      type = ErrorFractionMarker
      indicator = jump_von_mises
      refine = 0.7
      coarsen = 0.2
    []

    [boundary_marker]
      type = BoundaryMarker
      next_to = 'bottom'
      mark = refine
    []

    [boundary_marker_1]
      type = BoundaryMarker
      next_to = 'top'
      mark = refine
    []

    [combo]
      type = ComboMarker
      markers = 'errorfrac'
      # markers = 'boundary_marker boundary_marker_1'
      # markers = 'errorfrac boundary_marker boundary_marker_1'
    []
  []
[]

[keassyguang]

Hi @GiudGiud, sorry but do you have an update for the above non-working example?

[GiudGiud]

thanks for making the example. I'll be delayed but I won't forget.
If anyone else has time, please go ahead

[keassyguang]

Thank you very much, and I look forward to seeing the update.

[keassyguang]

I have written some debug output inside the InternalSideIndicator.C. What I found was:

1. the Real value = _field_var.dofValues()[0]; in the InternalSideIndicator::finalize() always returns zero when we set use_displaced_mesh = true as the global parameter.
2. I also checked the InternalSideIndicator::computeIndicator(), the for loop seems to be correct, however, the _field_var.nodalDofIndex() and _field_var.nodalDofIndexNeighbor() always point to 0. Not sure if it is reasonable.

I'm not sure whether it's the cause of the problem I've found before, also I don't know what to do next, any suggestions to fix this bug?

InternalSideIndicator.C

//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "InternalSideIndicator.h"

// MOOSE includes
#include "Assembly.h"
#include "MooseTypes.h"
#include "MooseVariableFE.h"
#include "Problem.h"
#include "SubProblem.h"
#include "SystemBase.h"

#include "libmesh/dof_map.h"
#include "libmesh/dense_vector.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_subvector.h"
#include "libmesh/libmesh_common.h"
#include "libmesh/quadrature.h"

const BoundaryID InternalSideIndicator::InternalBndId = 12345;

InputParameters
InternalSideIndicator::validParams()
{
  InputParameters params = Indicator::validParams();
  params.addRequiredParam<VariableName>(
      "variable", "The name of the variable that this side indicator applies to");
  params.addParam<bool>("scale_by_flux_faces",
                        false,
                        "Whether or not to scale the error values by "
                        "the number of flux faces.  This attempts to "
                        "not penalize elements on boundaries for "
                        "having less neighbors.");

  params.addPrivateParam<BoundaryID>("_boundary_id", InternalSideIndicator::InternalBndId);
  return params;
}

InternalSideIndicator::InternalSideIndicator(const InputParameters & parameters)
  : Indicator(parameters),
    NeighborCoupleable(this, false, false),
    ScalarCoupleable(this),
    NeighborMooseVariableInterface(
        this, false, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD),
    _field_var(_subproblem.getStandardVariable(_tid, name())),

    _current_elem(_assembly.elem()),
    _neighbor_elem(_assembly.neighbor()),

    _current_side(_assembly.side()),
    _current_side_elem(_assembly.sideElem()),

    _coord_sys(_assembly.coordSystem()),
    _q_point(_assembly.qPointsFace()),
    _qrule(_assembly.qRuleFace()),
    _JxW(_assembly.JxWFace()),
    _coord(_assembly.coordTransformation()),

    _boundary_id(parameters.get<BoundaryID>("_boundary_id")),

    _var(mooseVariableField()),
    _scale_by_flux_faces(parameters.get<bool>("scale_by_flux_faces")),

    _u(_var.sln()),
    _grad_u(_var.gradSln()),

    _normals(_assembly.normals()),

    _u_neighbor(_var.slnNeighbor()),
    _grad_u_neighbor(_var.gradSlnNeighbor())
{
  const std::vector<MooseVariableFieldBase *> & coupled_vars = getCoupledMooseVars();
  for (const auto & var : coupled_vars)
    addMooseVariableDependency(var);

  addMooseVariableDependency(&mooseVariableField());
}

void
InternalSideIndicator::computeIndicator()
{
  Real sum = 0;

  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
    sum += _JxW[_qp] * _coord[_qp] * computeQpIntegral();

  {
    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);

    // We use the nodal dof index but these are constant monomial variables
    _solution.add(_field_var.nodalDofIndex(), sum * _current_elem->hmax());
    _solution.add(_field_var.nodalDofIndexNeighbor(), sum * _neighbor_elem->hmax());
  }
}

void
InternalSideIndicator::finalize()
{
  unsigned int n_flux_faces = 0;

  if (_scale_by_flux_faces)
  {
    if (_var.isFV())
      paramError("scale_by_flux_faces", "Unsupported at this time for finite volume variables");

    // Figure out the total number of sides contributing to the error.
    // We'll scale by this so boundary elements are less penalized
    for (unsigned int side = 0; side < _current_elem->n_sides(); side++)
      if (_current_elem->neighbor_ptr(side) != nullptr)
        n_flux_faces++;
  }
  else
    n_flux_faces = 1;

  // The 0 is because CONSTANT MONOMIALS only have one coefficient per element
  Real value = _field_var.dofValues()[0];

  {
    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
    // We use the nodal dof index but these are constant monomial variables
    _solution.set(_field_var.nodalDofIndex(), std::sqrt(value) / static_cast<Real>(n_flux_faces));
  }
}