Isogeometric Analysis in MOOSE

This example illustrates using Isogeometric Analysis (IGA) within MOOSE framework to perform a simulation. This example simulates loads applied to a "c-frame" to determine the maximum principal stress.

Creating IGA Mesh

To create the mesh file included in this example a Coreform Cubit license is needed. Coreform Cubit is a product released and maintained by Coreform, LLC. A free to use Cubit Learn license can be acquired.

note:Mesh Generation with Cubit

The following is not required to run the example, but is required to generate the input mesh that is included in the example. To execute the Cubit journal file in batch from the command line use the following command:


coreform_cubit -batch cframe_build.jou

Within the journal file, see Listing 1, there are unique commands that will generate a uspline on the discretized mesh.

Line 35 sets the degree and continuity of the uspline
As of this writing, the max degree supported by libMesh is 2.
All continuity must equal p-1 where p is the degree.
Line 36 constructs the uspline using the geometry as a basis.
Line 37 fits the uspline that was built to the geometry.

Listing 1: Complete Coreform Cubit file for generating IGA input mesh

import step "c-frame.step" heal
remove surface 15 17 21 19 22 23 16 25 18 20 24 26 14 35 38 33 34 36 37 7 8 4 5 6 1 2 10 12 3 9 13 11 80 77 81 88 79 86 90 78 76 89 87 75 91 85 29 28 30 27 31 32 62 60 57 58 48 45 47 50 49 46 54 53 51 56 52 55 41 39 42 43 40 44 68 66  extend 
remove surface 93  extend 
remove surface 94  extend 
webcut volume 1  with sheet extended from surface 99 98 96  preview 
webcut volume 1  with sheet extended from surface 99 98 96  
webcut volume all with sheet extended from surface 92 59  
webcut volume 6 4 3 5 1 2  with sheet extended from surface 67 61  
webcut volume 1 16 2 15  with plane vertex 52  vertex 215  vertex 10  preview 
webcut volume 1 16 2 15  with plane vertex 52  vertex 215  vertex 10  
create vertex on curve 14  fraction 0.5 from start 
create vertex on curve 77  fraction 0.5 from start 
create vertex on curve 32  fraction 0.5 from start 
webcut volume 16 2 15  with plane vertex 270  vertex 269  vertex 271  
webcut volume 1  with plane vertex 270  vertex 269  vertex 271  preview 
webcut volume 1  with plane vertex 270  vertex 269  vertex 271  
delete Vertex 271 269 270
imprint all 
merge all
volume all redistribute nodes off 
volume all scheme Sweep Vector 1 0 0 
volume all scheme Sweep Vector 1 0 0 sweep transform least squares 
volume all autosmooth target on  fixed imprints off  smart smooth off 
volume all size 0.5
mesh volume all
block 1 vol all
sideset 2 surf 194 203 196 
sideset 3 Surface 146 155 148 
block 1 element type hex27 #Quadratic FEA elements
export mesh "CFrame_FEA_05_deg2.e"  overwrite #FEA Export
set uspline vol all deg 2 cont 1
build uspline vol all as 1
fit uspline 1
export uspline 1 exodus "CFrame_IGA_05"

MOOSE-IGA Simulation

Performing the simulation utilizing the mesh created above does not require much with respect to the MOOSE input, simply load the mesh from a file and select utilize the RATIONAL_BERNSTEIN element family as shown in Listing 2. Exporting using the VTK format (vtk = true) input will output in a format that will capture the higher-order nature of the IGA based elements using Paraview visualization.

Listing 2: Complete input file for running example problem with IGA in MOOSE.

[GlobalParams<<<{"href": "../../../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y disp_z'
[]

[Mesh<<<{"href": "../../../syntax/Mesh/index.html"}>>>]
  [igafile]
    type = FileMeshGenerator<<<{"description": "Read a mesh from a file.", "href": "../../../source/meshgenerators/FileMeshGenerator.html"}>>>
    file<<<{"description": "The filename to read."}>>> = cframe_iga_coarse.e
    clear_spline_nodes<<<{"description": "If clear_spline_nodes=true, IsoGeometric Analyis spline nodes and constraints are removed from an IGA mesh, after which only C^0 Rational-Bernstein-Bezier elements will remain."}>>> = true
  []
[]

[Variables<<<{"href": "../../../syntax/Variables/index.html"}>>>]
  [disp_x]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = SECOND
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = RATIONAL_BERNSTEIN
  []
  [disp_y]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = SECOND
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = RATIONAL_BERNSTEIN
  []
  [disp_z]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = SECOND
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = RATIONAL_BERNSTEIN
  []
[]

[Kernels<<<{"href": "../../../syntax/Kernels/index.html"}>>>]
  [SolidMechanics<<<{"href": "../../../syntax/Kernels/SolidMechanics/index.html"}>>>]
#Stress divergence kernels
    displacements<<<{"description": "The nonlinear displacement variables for the problem"}>>> = 'disp_x disp_y disp_z'
   []
[]

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
    [von_mises]
        #Dependent variable used to visualize the von Mises stress
       order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = SECOND
       family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    []
    [Max_Princ]
       #Dependent variable used to visualize the Hoop stress
       order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = SECOND
       family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    []
    [stress_xx]
        order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = SECOND
        family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    []
    [stress_yy]
        order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = SECOND
        family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    []
    [stress_zz]
        order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = SECOND
        family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    []
[]

[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [von_mises_kernel]
    #Calculates the von mises stress and assigns it to von_mises
    type = RankTwoScalarAux<<<{"description": "Compute a scalar property of a RankTwoTensor", "href": "../../../source/auxkernels/RankTwoScalarAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = von_mises
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    scalar_type<<<{"description": "Type of scalar output"}>>> = VonMisesStress
  []
  [MaxPrin]
    type = RankTwoScalarAux<<<{"description": "Compute a scalar property of a RankTwoTensor", "href": "../../../source/auxkernels/RankTwoScalarAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = Max_Princ
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    scalar_type<<<{"description": "Type of scalar output"}>>> = MaxPrincipal
  []
  [stress_xx]
    type = RankTwoAux<<<{"description": "Access a component of a RankTwoTensor", "href": "../../../source/auxkernels/RankTwoAux.html"}>>>
    index_i<<<{"description": "The index i of ij for the tensor to output (0, 1, 2)"}>>> = 0
    index_j<<<{"description": "The index j of ij for the tensor to output (0, 1, 2)"}>>> = 0
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_xx
  []
    [stress_yy]
    type = RankTwoAux<<<{"description": "Access a component of a RankTwoTensor", "href": "../../../source/auxkernels/RankTwoAux.html"}>>>
    index_i<<<{"description": "The index i of ij for the tensor to output (0, 1, 2)"}>>> = 1
    index_j<<<{"description": "The index j of ij for the tensor to output (0, 1, 2)"}>>> = 1
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_yy
  []
  [stress_zz]
    type = RankTwoAux<<<{"description": "Access a component of a RankTwoTensor", "href": "../../../source/auxkernels/RankTwoAux.html"}>>>
    index_i<<<{"description": "The index i of ij for the tensor to output (0, 1, 2)"}>>> = 2
    index_j<<<{"description": "The index j of ij for the tensor to output (0, 1, 2)"}>>> = 2
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_zz
  []
[]

[BCs<<<{"href": "../../../syntax/BCs/index.html"}>>>]
  [Pressure<<<{"href": "../../../syntax/BCs/Pressure/index.html"}>>>]
    [load]
      #Applies the pressure
      boundary<<<{"description": "The list of boundary IDs from the mesh where the pressure will be applied"}>>> = '3'
      factor<<<{"description": "The factor to use in computing the pressure"}>>> = 2000 # psi
    []
  []
  [anchor_x]
    #Anchors the bottom and sides against deformation in the x-direction
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = '2'
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [anchor_y]
    #Anchors the bottom and sides against deformation in the y-direction
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = '2'
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [anchor_z]
    #Anchors the bottom and sides against deformation in the z-direction
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_z
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = '2'
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
[]

[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
  [elasticity_tensor_AL]
    #Creates the elasticity tensor using concrete parameters
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 24e6 #psi
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0.33
    type = ComputeIsotropicElasticityTensor<<<{"description": "Compute a constant isotropic elasticity tensor.", "href": "../../../source/materials/ComputeIsotropicElasticityTensor.html"}>>>
  []
  [strain]
    #Computes the strain, assuming small strains
    type = ComputeSmallStrain<<<{"description": "Compute a small strain.", "href": "../../../source/materials/ComputeSmallStrain.html"}>>>
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
  [stress]
    #Computes the stress, using linear elasticity
    type = ComputeLinearElasticStress<<<{"description": "Compute stress using elasticity for small strains", "href": "../../../source/materials/ComputeLinearElasticStress.html"}>>>
  []
  [density_AL]
    #Defines the density of steel
    type = GenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../../../source/materials/GenericConstantMaterial.html"}>>>
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = density
    prop_values<<<{"description": "The values associated with the named properties"}>>> = 6.99e-4 # lbm/in^3
  []
[]

[Preconditioning<<<{"href": "../../../syntax/Preconditioning/index.html"}>>>]
  [SMP]
    #Creates the entire Jacobian, for the Newton solve
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
  []
[]

[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [max_principal_stress]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0.000000 -1.500000 -4.3'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = Max_Princ
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation.  Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
  []
  [maxPrincStress]
    type = ElementExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../../../source/postprocessors/ElementExtremeValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = Max_Princ
  []
[]

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  # We solve a steady state problem using Newton's iteration
  type = Steady
  solve_type = NEWTON
  nl_rel_tol = 1e-9
  l_max_its = 300
  l_tol = 1e-4
  nl_max_its = 30
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 31'
[]

[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
  vtk<<<{"description": "Output the results using the default settings for VTKOutput output"}>>> = true
[]

Figure 1: Maximum principal stress for "c-frame" example utilizing IGA in MOOSE.

(moose/modules/solid_mechanics/examples/cframe_iga/cframe_build.jou)

import step "c-frame.step" heal
remove surface 15 17 21 19 22 23 16 25 18 20 24 26 14 35 38 33 34 36 37 7 8 4 5 6 1 2 10 12 3 9 13 11 80 77 81 88 79 86 90 78 76 89 87 75 91 85 29 28 30 27 31 32 62 60 57 58 48 45 47 50 49 46 54 53 51 56 52 55 41 39 42 43 40 44 68 66  extend 
remove surface 93  extend 
remove surface 94  extend 
webcut volume 1  with sheet extended from surface 99 98 96  preview 
webcut volume 1  with sheet extended from surface 99 98 96  
webcut volume all with sheet extended from surface 92 59  
webcut volume 6 4 3 5 1 2  with sheet extended from surface 67 61  
webcut volume 1 16 2 15  with plane vertex 52  vertex 215  vertex 10  preview 
webcut volume 1 16 2 15  with plane vertex 52  vertex 215  vertex 10  
create vertex on curve 14  fraction 0.5 from start 
create vertex on curve 77  fraction 0.5 from start 
create vertex on curve 32  fraction 0.5 from start 
webcut volume 16 2 15  with plane vertex 270  vertex 269  vertex 271  
webcut volume 1  with plane vertex 270  vertex 269  vertex 271  preview 
webcut volume 1  with plane vertex 270  vertex 269  vertex 271  
delete Vertex 271 269 270
imprint all 
merge all
volume all redistribute nodes off 
volume all scheme Sweep Vector 1 0 0 
volume all scheme Sweep Vector 1 0 0 sweep transform least squares 
volume all autosmooth target on  fixed imprints off  smart smooth off 
volume all size 0.5
mesh volume all
block 1 vol all
sideset 2 surf 194 203 196 
sideset 3 Surface 146 155 148 
block 1 element type hex27 #Quadratic FEA elements
export mesh "CFrame_FEA_05_deg2.e"  overwrite #FEA Export
set uspline vol all deg 2 cont 1
build uspline vol all as 1
fit uspline 1
export uspline 1 exodus "CFrame_IGA_05"

(moose/modules/solid_mechanics/examples/cframe_iga/cframe_iga.i)

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Mesh]
  [igafile]
    type = FileMeshGenerator
    file = cframe_iga_coarse.e
    clear_spline_nodes = true
  []
[]

[Variables]
  [disp_x]
    order = SECOND
    family = RATIONAL_BERNSTEIN
  []
  [disp_y]
    order = SECOND
    family = RATIONAL_BERNSTEIN
  []
  [disp_z]
    order = SECOND
    family = RATIONAL_BERNSTEIN
  []
[]

[Kernels]
  [SolidMechanics]
#Stress divergence kernels
    displacements = 'disp_x disp_y disp_z'
   []
[]

[AuxVariables]
    [von_mises]
        #Dependent variable used to visualize the von Mises stress
       order = SECOND
       family = MONOMIAL
    []
    [Max_Princ]
       #Dependent variable used to visualize the Hoop stress
       order = SECOND
       family = MONOMIAL
    []
    [stress_xx]
        order = SECOND
        family = MONOMIAL
    []
    [stress_yy]
        order = SECOND
        family = MONOMIAL
    []
    [stress_zz]
        order = SECOND
        family = MONOMIAL
    []
[]

[AuxKernels]
  [von_mises_kernel]
    #Calculates the von mises stress and assigns it to von_mises
    type = RankTwoScalarAux
    variable = von_mises
    rank_two_tensor = stress
    scalar_type = VonMisesStress
  []
  [MaxPrin]
    type = RankTwoScalarAux
    variable = Max_Princ
    rank_two_tensor = stress
    scalar_type = MaxPrincipal
  []
  [stress_xx]
    type = RankTwoAux
    index_i = 0
    index_j = 0
    rank_two_tensor = stress
    variable = stress_xx
  []
    [stress_yy]
    type = RankTwoAux
    index_i = 1
    index_j = 1
    rank_two_tensor = stress
    variable = stress_yy
  []
  [stress_zz]
    type = RankTwoAux
    index_i = 2
    index_j = 2
    rank_two_tensor = stress
    variable = stress_zz
  []
[]

[BCs]
  [Pressure]
    [load]
      #Applies the pressure
      boundary = '3'
      factor = 2000 # psi
    []
  []
  [anchor_x]
    #Anchors the bottom and sides against deformation in the x-direction
    type = DirichletBC
    variable = disp_x
    boundary = '2'
    value = 0.0
  []
  [anchor_y]
    #Anchors the bottom and sides against deformation in the y-direction
    type = DirichletBC
    variable = disp_y
    boundary = '2'
    value = 0.0
  []
  [anchor_z]
    #Anchors the bottom and sides against deformation in the z-direction
    type = DirichletBC
    variable = disp_z
    boundary = '2'
    value = 0.0
  []
[]

[Materials]
  [elasticity_tensor_AL]
    #Creates the elasticity tensor using concrete parameters
    youngs_modulus = 24e6 #psi
    poissons_ratio = 0.33
    type = ComputeIsotropicElasticityTensor
  []
  [strain]
    #Computes the strain, assuming small strains
    type = ComputeSmallStrain
    displacements = 'disp_x disp_y disp_z'
  []
  [stress]
    #Computes the stress, using linear elasticity
    type = ComputeLinearElasticStress
  []
  [density_AL]
    #Defines the density of steel
    type = GenericConstantMaterial
    prop_names = density
    prop_values = 6.99e-4 # lbm/in^3
  []
[]

[Preconditioning]
  [SMP]
    #Creates the entire Jacobian, for the Newton solve
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [max_principal_stress]
    type = PointValue
    point = '0.000000 -1.500000 -4.3'
    variable = Max_Princ
    use_displaced_mesh = false
  []
  [maxPrincStress]
    type = ElementExtremeValue
    variable = Max_Princ
  []
[]

[Executioner]
  # We solve a steady state problem using Newton's iteration
  type = Steady
  solve_type = NEWTON
  nl_rel_tol = 1e-9
  l_max_its = 300
  l_tol = 1e-4
  nl_max_its = 30
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 31'
[]

[Outputs]
  vtk = true
[]