adaptive mesh refinement and block generation

[keassyguang]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Hello everyone,
Recently, I found a paper dealing with the additive manufacturing process simulation with adaptive mesh refinement implemented via DEAL.II. The key adaptive mesh refinement concept in that paper is to use a coarse mesh for the initial set up and then refine/coarsen the mesh as needed to reduce computational cost. As shown in the following picture, 1) the top-most active layer contains 1 element in the build direction; 2) elements near the top surface are refined; 3) elements far away from the top surface are coarsen; 4) elements near the part boundaries are refined. I want to implement this concept with MOOSE. However, I got some problems.
Suppose a block of size 4e-3 x 4e-3 x 4e-3 is initially meshed with 5x5x5 elements. The layer thickness is 50e-6, so that 4 refinement levels are needed to reach such a cell size. The top surface moves every second as defined by CoupledVarThresholdElementSubdomainModifier and an AuxVariable. An error message is displayed: "The block '1' has not been found in the network" when I run the code. I know that if I use 80 (4e-3 / 50e-6 = 80) elements, the code can run. However, the initial mesh in such a case is already a fine mesh and I can't coarsen the mesh to a larger size. How can I get the mesh to look like the picture?

[Problem]
  kernel_coverage_check = false
  boundary_restricted_node_integrity_check = false
  boundary_restricted_elem_integrity_check = false
  type = FEProblem
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 5
    ny = 5
    nz = 5
    xmax = 4e-3
    ymax = 4e-3
    zmax = 4e-3
  []
  [active_domain] #the initially active domain
    type = SubdomainBoundingBoxGenerator
    input = 'gen'
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '4e-3 4e-3 50e-6'
  []
  [inactive_domain] #the initially inactive domain
    type = SubdomainBoundingBoxGenerator
    input = active_domain
    block_id = 2
    bottom_left = '0 0 50e-6'
    top_right = '4e-3 4e-3 4e-3'
  []
  [sideset_top]
    input = inactive_domain
    type = SideSetsAroundSubdomainGenerator
    normal = '0 0 1'
    block = 1
    new_boundary = 'sideset_top' # top boundary of active domain
  []
[]

[Variables]
  [T]
    initial_condition = 300 #Kelvin
    block = 1
  []
[]

[AuxVariables]
  [layerZ]
    [AuxKernel]
      type = ParsedAux
      expression = 'z - 50e-6 * ceil(t/1.000001)'
      use_xyzt = true
      execute_on = 'INITIAL TIMESTEP_BEGIN'
    []
  []
[] 

[Kernels]
  [conduction]
    block = 1
    type = ADHeatConduction
    variable = T
  []
[]

[MeshModifiers]
  [addLayer] 
    type = CoupledVarThresholdElementSubdomainModifier
    coupled_var = 'layerZ'
    criterion_type = BELOW
    threshold = 0
    subdomain_id = 1
    moving_boundaries = 'sideset_top'
    moving_boundary_subdomain_pairs = '1 2'
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
[]

[Materials]
  [thermal]
    block = 1
    type = ADHeatConductionMaterial
    specific_heat = 1
    thermal_conductivity = 1e-3
  []
  [density]
    block = 1
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = '1'
  []
  [dummy_mat_inactive]
    block = 2
    type = GenericConstantMaterial
    prop_names = 'dummy_mat'
    prop_values = '0.0'
  []
[]

[BCs]
  # In 3D, back = 0, bottom = 1, right = 2, top = 3, left = 4, front = 5
  [bottom]
    type = ADDirichletBC
    variable = T
    boundary = 0
    value = 400
  []
  [side_top_BC]
    type = ADConvectiveHeatFluxBC
    T_infinity = 500
    heat_transfer_coefficient = 100
    variable = T
    boundary = sideset_top
  []
[]

[Executioner]
  type = Transient
  solve_type = PJFNK
  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type -pc_hypre_boomeramg_max_iter'
  petsc_options_value = '201                hypre    boomeramg      8'
  end_time = 80
  dt = 0.1
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-12
[]

[Outputs]
  exodus = true
[]

[GiudGiud]

Does your starting mesh show block 1 as expected?

Which object is returning the error? It should say in the console

This combination of adaptivity criteria should be easy to do in moose. We can combine Markers to target each of these needs

[keassyguang]

I simplified the input file to a 2D case, and here is the initial mesh generated within below input file. However, the active and inactive domains and the moving boundary are not as expected as shown in the picture. Also, the moving boundary moves in an incorrect way. I would like the moving boundary to move along the Y-direction 50e-6 every second.

1. What should I do to make the active and inactive domains and the moving boundary change correctly?
2. Which marker is suitable for refining or coarsening the mesh near a boundary? Here I want to refine the mesh near the moving boundary and coarsen the mesh far from the moving boundary.

[Problem]
  kernel_coverage_check = false
  boundary_restricted_node_integrity_check = false
  boundary_restricted_elem_integrity_check = false
  type = FEProblem
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 5
    ny = 5
    xmax = 4e-3
    ymax = 4e-3
  []
  [active_domain] #the initially active domain
    type = SubdomainBoundingBoxGenerator
    input = 'gen'
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '4e-3 800e-6 0'
    # top_right = '4e-3 50e-6 0'
  []
  [inactive_domain] #the initially inactive domain
    type = SubdomainBoundingBoxGenerator
    input = active_domain
    block_id = 2
    bottom_left = '0 800e-6 0'
    # bottom_left = '0 50e-6 0'
    top_right = '4e-3 4e-3 0'
  []
  # [sideset_top]
  #   input = inactive_domain
  #   type = SideSetsAroundSubdomainGenerator
  #   normal = '0 1 0'
  #   block = 1
  #   new_boundary = 'sideset_top' # top boundary of active domain
  # []
  [sideset_top]
    type = SideSetsBetweenSubdomainsGenerator
    input = inactive_domain
    new_boundary = 'sideset_top'
    paired_block = 2
    primary_block = 1
  []
[]

[Variables]
  [T]
    initial_condition = 300 #Kelvin
    block = 1
  []
[]

[AuxVariables]
  [layerY]
    [AuxKernel]
      type = ParsedAux
      expression = 'y - 50e-6 * ceil(t/1.000001)'
      use_xyzt = true
      execute_on = 'INITIAL TIMESTEP_BEGIN'
    []
  []
[] 

[Kernels]
  [conduction]
    block = 1
    type = ADHeatConduction
    variable = T
  []
[]

[MeshModifiers]
  [addLayer] 
    type = CoupledVarThresholdElementSubdomainModifier
    coupled_var = 'layerY'
    criterion_type = BELOW
    threshold = 0
    subdomain_id = 1
    moving_boundaries = 'sideset_top'
    moving_boundary_subdomain_pairs = '1 2'
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
[]

[Materials]
  [thermal]
    block = 1
    type = ADHeatConductionMaterial
    specific_heat = 1
    thermal_conductivity = 1e-3
  []
  [density]
    block = 1
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = '1'
  []
  [dummy_mat_inactive]
    block = 2
    type = GenericConstantMaterial
    prop_names = 'dummy_mat'
    prop_values = '0.0'
  []
[]

[BCs]
  [bottom]
    type = ADDirichletBC
    variable = T
    boundary = 0
    value = 400
  []
  [side_top_BC]
    type = ADConvectiveHeatFluxBC
    T_infinity = 500
    heat_transfer_coefficient = 100
    variable = T
    boundary = sideset_top
  []
[]

[Executioner]
  type = Transient
  solve_type = PJFNK
  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type -pc_hypre_boomeramg_max_iter'
  petsc_options_value = '201                hypre    boomeramg      8'
  end_time = 80
  dt = 0.1
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-12
[]

[Outputs]
  exodus = true
[]

[Adaptivity]
  [Markers]
    [boundary]
      type = BoundaryMarker
      next_to = 'bottom'
      mark = refine
    []
    [boundary1]
      type = BoundaryMarker
      next_to = 'sideset_top'
      mark = refine
    []
    # [box]
    #   type = BoxMarker
    #   bottom_left = '0 y-50e-6*ceil(t/1.000001)+50e-6 0'
    #   top_right = '4e-3 y-50e-6*ceil(t/1.000001)-3*50e-6 0'
    #   inside = refine
    #   outside = coarsen
    # []
  []
  initial_marker = boundary
  initial_steps = 4
  max_h_level = 4
  marker = boundary1
[]

[GiudGiud]

What should I do to make the active and inactive domains and the moving boundary change correctly?

For moving subdomains and boundaries, you can use this object
https://mooseframework.inl.gov/source/meshmodifiers/CoupledVarThresholdElementSubdomainModifier.html
then for example, using a FunctionAux for example to make an auxiliary variable describe the subdomain changes to make happen

Which marker is suitable for refining or coarsening the mesh near a boundary? Here I want to refine the mesh near the moving boundary and coarsen the mesh far from the moving boundary.

This one
https://mooseframework.inl.gov/source/markers/BoundaryMarker.html

[keassyguang]

Thank you very much. The objects I used are exactly the same as the ones you mentioned. However, the behaviour of them is not correct. Can you kindly help me check the input file to find out where I am going wrong?

[Problem]
  kernel_coverage_check = false
  boundary_restricted_node_integrity_check = false
  boundary_restricted_elem_integrity_check = false
  type = FEProblem
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 5
    ny = 5
    xmax = 4e-3
    ymax = 4e-3
  []
  [active_domain] #the initially active domain
    type = SubdomainBoundingBoxGenerator
    input = 'gen'
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '4e-3 800e-6 0'
    # top_right = '4e-3 50e-6 0'
  []
  [inactive_domain] #the initially inactive domain
    type = SubdomainBoundingBoxGenerator
    input = active_domain
    block_id = 2
    bottom_left = '0 800e-6 0'
    # bottom_left = '0 50e-6 0'
    top_right = '4e-3 4e-3 0'
  []
  # [sideset_top]
  #   input = inactive_domain
  #   type = SideSetsAroundSubdomainGenerator
  #   normal = '0 1 0'
  #   block = 1
  #   new_boundary = 'sideset_top' # top boundary of active domain
  # []
  [sideset_top]
    type = SideSetsBetweenSubdomainsGenerator
    input = inactive_domain
    new_boundary = 'sideset_top'
    paired_block = 2
    primary_block = 1
  []
[]

[Variables]
  [T]
    initial_condition = 300 #Kelvin
    block = 1
  []
[]

[AuxVariables]
  [layerY]
    [AuxKernel]
      type = ParsedAux
      expression = 'y - 50e-6 * ceil(t/1.000001)'
      use_xyzt = true
      execute_on = 'INITIAL TIMESTEP_BEGIN'
    []
  []
[] 

[Kernels]
  [conduction]
    block = 1
    type = ADHeatConduction
    variable = T
  []
[]

[MeshModifiers]
  [addLayer] 
    type = CoupledVarThresholdElementSubdomainModifier
    coupled_var = 'layerY'
    criterion_type = BELOW
    threshold = 0
    subdomain_id = 1
    moving_boundaries = 'sideset_top'
    moving_boundary_subdomain_pairs = '1 2'
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
[]

[Materials]
  [thermal]
    block = 1
    type = ADHeatConductionMaterial
    specific_heat = 1
    thermal_conductivity = 1e-3
  []
  [density]
    block = 1
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = '1'
  []
  [dummy_mat_inactive]
    block = 2
    type = GenericConstantMaterial
    prop_names = 'dummy_mat'
    prop_values = '0.0'
  []
[]

[BCs]
  [bottom]
    type = ADDirichletBC
    variable = T
    boundary = 0
    value = 400
  []
  [side_top_BC]
    type = ADConvectiveHeatFluxBC
    T_infinity = 500
    heat_transfer_coefficient = 100
    variable = T
    boundary = sideset_top
  []
[]

[Executioner]
  type = Transient
  solve_type = PJFNK
  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type -pc_hypre_boomeramg_max_iter'
  petsc_options_value = '201                hypre    boomeramg      8'
  end_time = 80
  dt = 0.1
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-12
[]

[Outputs]
  exodus = true
[]

[Adaptivity]
  [Markers]
    [boundary]
      type = BoundaryMarker
      next_to = 'bottom'
      mark = refine
    []
    [boundary1]
      type = BoundaryMarker
      next_to = 'sideset_top'
      mark = refine
    []
    # [box]
    #   type = BoxMarker
    #   bottom_left = '0 y-50e-6*ceil(t/1.000001)+50e-6 0'
    #   top_right = '4e-3 y-50e-6*ceil(t/1.000001)-3*50e-6 0'
    #   inside = refine
    #   outside = coarsen
    # []
  []
  initial_marker = boundary
  initial_steps = 4
  max_h_level = 4
  marker = boundary1
[]

[GiudGiud]

I can take a look tomorrow.
What part is not working as you expect? The boundary movement or the refinement?

I think the coarsening once the boundary has moved wont be automatic. Or if it is, it will have to be handled by a separate Marker not by the boundary marker

[keassyguang]

Thanks a lot. Both the boundary movement and refinement do not work as expected. Sometimes they both work correctly, but sometimes they do not. Also, I can't understand why so many boundaries are generated.

I think the coarsening once the boundary has moved wont be automatic. Or if it is, it will have to be handled by a separate Marker not by the boundary marker

Correct. I need additional markers to coarsen the mesh far away from the moving boundary.
Below picture shows some results based on the above input file.

[GiudGiud]

Both the boundary movement

It's on a per element basis, and the criteria is a comparison to 0 which is a value that the criteria hits exactly. So that's why you see dents like at 26.1

ceil(t/1.000001)

is not a great idea in the expression. This means that sometimes you get two (or more) steps that have the same value for ceil
so then the criteria does not move even though the time changed

[keassyguang]

I tried execute_on = 'TIMESTEP_END', but it seems the results are the same. I also tried another marker, ValueRangeMarker. The bottom most elements of block 2 (red color) were designed to refine to max_h_level, then the mesh size of the top elements of block 1 would always be at the max_h_level when the moving boundary moves to a new position.

As for the second problem, I didn't see any possible solution. It seems that it cannot be achieved by MOOSE and DEALII either.

[GiudGiud]

The results are not the same. If you look at what you posted a week ago, the most refined area is now right near the boundary, when it was 1 refinement step removed.

Now if you want to have the red area be most refined, you need to use the sideset oriented the other direction. Right now, it is refining below at the highest level. The sideset oriented the other way will refine the red first.

As for the second problem, I didn't see any possible solution. It seems that it cannot be achieved by MOOSE and DEALII either.

I explained a solution, you just need to code a custom Marker. MOOSE lets you code the exact behavior you want. If you want no-code solutions you ll have to find an app that has been developed to study your problem

[keassyguang]

The results are not the same. If you look at what you posted a week ago, the most refined area is now right near the boundary, when it was 1 refinement step removed.

Yes, I want the most refined area is near the boundary. The results so far seem to meet my needs.

Now if you want to have the red area be most refined, you need to use the sideset oriented the other direction. Right now, it is refining below at the highest level. The sideset oriented the other way will refine the red first.

Thanks for your explanation and suggestions.

You could try to prevent the coarsening cells that "contain a boundary" using a custom marker.
I explained a solution, you just need to code a custom Marker.

Do you mean this solution? Since the top elements of block 1 (white colour) have already reached the highest level of refinement, if we prevent the elements containing the boundary from being coarsened, these elements will remain the highest refinement until the boundary moves to the next position. What I want is to refine elements near the boundary but inside the manufactured part. In this case I have no idea what to do.

[GiudGiud]

You need to replace the combination of two markers that does this

Since the top elements of block 1 (white colour) have already reached the highest level of refinement, if we prevent the elements containing the boundary from being coarsened, these elements will remain the highest refinement until the boundary moves to the next position.

with a custom marker that does this

What I want is to refine elements near the boundary but inside the manufactured part

Markers are fairly simple objects to code, same as kernels or auxkernels. You have a function to implement that dictates the refine/coarsen decision on an element

[keassyguang]

Thank you very much. I may explore it a bit later as I need to examine the mechanical model to make sure the whole framework works.

BTW, do we have a user object to output the number of elements of the model at each time step?