# 14.2. Usage Guidelines

The following block command defines a `CRK SECTION`.

    # Within Sierra scope...
    BEGIN CRK SECTION <string>crk_name
      # Commands controlling approximation properties
      SHAPE FUNCTION TYPE = SPHERICAL|SPHERICAL_SURFACE_KRONECKER|
                            TWO_STAR|LINEAR_TET (SPHERICAL)
      SUPPORT SIZE = <real>support_size (1.1)
      SUPPORT SIZE TYPE = RELATIVE_LOCAL|RELATIVE_GLOBAL|
                          PHYSICAL (RELATIVE_LOCAL)
      BLOCK CONFORMING = ON|OFF (ON)

      # Commands controlling numerical integration
      FORMULATION = ELEMENT|NODAL_FROM_THEX|TET_FROM_THEX (ELEMENT)
      NUMBER OF INTEGRATION SUBCELLS = <int>num_sub_cells (1)

      # Controls which stress-strain conjugate pairs to use
      # ADAGIO uses a Cauchy-stress/rate-of-deformation conjugate
      # LOGARITHMIC uses Hencky strain/"Hencky" stress
      STRAIN MEASURE = ADAGIO|LOGARITHMIC (ADAGIO)

      # Controls when post processing occurs for performance
      POST PROCESS EVERY STEP = FALSE|TRUE (TRUE)

      # Creates face entities to apply contact to.  Required when
      # using Dash contact with CRK.
      CREATE CONTACT FACES = ON|OFF (OFF)

      # Commands controlling bond-based failure modeling in CRK
      FAILURE METHOD = BOND|NONE (NONE)
      FAILURE INDEX = <integer>failure_index
      FAILURE CRITERION = <real>bond_failure_criterion
      BOND OPENING DISPLACEMENT = <real>bond_opening_displacement
      BOND SECONDARY CRITERION = <real>bond_secondary_criterion
      SHOW BONDS = ON|OFF (OFF)

      # Controls which nodes are DOFs in the model
      # (for nearly-incompressible media)
      THEX ACTIVE NODES = NODES|EDGES|FACES|CENTROIDS (NODES)

    END [CRK SECTION <string>crk_name]

The `SHAPE FUNCTION TYPE`, `SUPPORT SIZE` `SUPPORT SIZE TYPE`, and `BLOCK CONFORMING` options control the approximation properties of the CRK formulation.

The command `SHAPE_FUNCTION_TYPE` determines the kernel shape used in the CRK analysis. The option `SPHERICAL` is for the spherical-star where the elements that are partially or fully contained within a ball are included in the support. Similarly, the `SPHERICAL_SURFACE_KRONECKER` uses spherical-star kernels for the interior of the domain but with the support restricted to a one-star at surfaces, giving the same shape functions as a linear tetrahedral finite element on the surface and thus possessing the Kronecker delta property at these locations. The option `TWO_STAR` uses kernels with a two-star kernel definition. Finally, the option `LINEAR_TET`, gives linear tetrahedral finite element shape functions, equivalent to using a one-star kernel on a tetrahedral mesh.

Warning

The option `SHAPE FUNCTION TYPE = TWO_STAR` is only compatible with `FORMULATION = ELEMENT` described below. Additionally, the input mesh must be a non `thex`, tetrahedral topology.

The commands `SUPPORT SIZE` and `SUPPORT SIZE TYPE` control the support size and are only valid for `SHAPE FUNCTION TYPE = SPHERICAL` or `SPHERICAL_SURFACE_KRONECKER`. `SUPPORT SIZE TYPE` dictates whether the size is relative to the maximum length of attached edges (`RELATIVE_LOCAL`), relative to the global max element edge length (`RELATIVE_GLOBAL`), or absolute (`PHYSICAL`). The command `SUPPORT SIZE` sets the absolute or relative support size. Relative support sizes near but greater than 1 are recommended.

The command `BLOCK CONFORMING = ON|OFF` provides the option to make the kernel functions to conform to block boundaries or not. Non-conforming will give a behavior similar to traditional meshfree methods.

The command `FORMULATION` determines which type of numerical integration should be performed and supports three options. The option `ELEMENT` integrates using the elements of the input mesh as the strain-smoothing domains. `NODAL_FROM_THEX` specifies that nodal integration/strain-smoothing domains should be used and requires tetrahedral elements that are subdivided beforehand into hexahedral elements using the `thex` command in Cubit. The nodal integration domains are formed by aggregating all of the sub-tet hex elements connected to each node to form the barycentric dual, referred to as an integration “cell.” This results in the SCNI strategy [[1]]. The option `TET_FROM_THEX` integrates over the original tetrahedra before subdivision. For the latter two options, the sub-tet hex elements can be used for sub-cells in the \(\mathbf{\overline{F}}\) integration technique as described in Section 14.2.2.

If `NODAL_FROM_THEX` or `TET_FROM_THEX` are specified above, the nodes in the CRK domain will be divided into an “active” group and an “inactive” group. The “active” group will be the DOF-carrying nodes and usually represents the original tet nodes before the `thex` command is issued in Cubit. The “inactive” group will contain all other nodes which are effectively evaluation points used to setup the nodal integration domains used in CRK. These complementary sets may be adjusted using the `THEX ACTIVE NODES` command as described in Section 14.2.2.

Given its nonlocal nature, the CRK method typically lacks the Kronecker delta property and thus incurs additional post-processing steps for output purposes. The `POST PROCESS EVERY STEP` command controls whether post-processing procedures occur every step (`TRUE`) or during output steps only (`FALSE`). This improves overall performance as these computations are not necessary for the internal force routine and will decrease the total runtime if set to `FALSE` (contingent on the output frequency). However, derived output or user routines which depend on CRK variables on every step will not function correctly as a result. As such, the command `POST PROCESS EVERY STEP = FALSE` should only be used when user output, user functions, and other derived variables do not depend on CRK fields. The default is `TRUE`.

The command `STRAIN MEASURE = ADAGIO|LOGARITHMIC` controls which strain to use when evaluating the constitutive law. The default is `ADAGIO` which uses the rate of deformation (the default in Thecode). The option `LOGARITHMIC` uses the logarithmic strain, or Hencky strain, and its corresponding conjugate stress instead.

## 14.2.1. Contact

CRK currently supports Dash contact but the setup will be slightly different depending on whether the command `FORMULATION` is set to `ELEMENT` versus `NODAL_FROM_THEX` or `TET_FROM_THEX`. For both cases, the shape function type must be either `SPHERICAL_SURFACE_KRONECKER` or `LINEAR_TET`.

In the case that `FORMULATION = ELEMENT`, every node in the mesh will be labeled as an “active” node and Dash contact may be setup as if the simulation was using normal elements without any additional modification as given in the Sierra/SM User Manual.

In the case of `FORMULATION = NODAL_FROM_THEX` or `TET_FROM_THEX`, the CRK surface must be manually created using the command `CREATE CONTACT FACES = ON`. This command will skin the CRK block(s) and create a new block with the name `<block_name>_contact_surface` where `<block_name>` is the original block name. This new surface is necessary to setup CRK shape function weights and area densities that Dash requires. These weights spread the contact forces from inactive nodes to active ones associated with the problem unknowns. Once this surface is created, it may be referenced within the contact definition as, for example

    # Within Sierra scope...
    BEGIN CRK SECTION crk_section
      SHAPE FUNCTION TYPE = SPHERICAL_SURFACE_KRONECKER
      FORMULATION = NODAL_FROM_THEX
      CREATE CONTACT FACES = ON
    END CRK SECTION crk_section
    ...
    # Within Region scope...
    BEGIN CONTACT DEFINITION collide
      SEARCH = Dash
      CONTACT SURFACE crk_surface CONTAINS block_1_contact_surface
      CONTACT SURFACE block_2 CONTAINS block_2

      BEGIN INTERACTION DEFAULTS
        GENERAL CONTACT = ON
        FRICTION MODEL = FRICTIONLESS
        INTERACTION BEHAVIOR = SLIDING
      END INTERACTION DEFAULTS
    END CONTACT DEFINITION collide

where in the above `block_1` is the CRK block and `block_2` is an FEM block. Dash will then take the information present on the CRK surface including shape function values and weights in order to create the contact representation.

## 14.2.2. Advanced integration methods

By default, the \(\mathbf{\overline{F}}\) integration technique [[2]] is activated if the command `FORMULATION =` `NODAL_FROM_THEX` or `TET_FROM_THEX` and if the `NUMBER OF INTEGRATION SUBCELLS` is greater than one. In this case, the sub-tet hex elements will be used for sub-cells in the \(\mathbf{\overline{F}}\) method for improving performance with nearly incompressible materials. In this context a “sub-cell” is a collection of sub-tet hex elements which share an original tet node and represents a deviatoric domain of the integration cell. The number of sub-cells is specified with `NUMBER OF INTEGRATION CELLS`. Elements will be grouped into sub-cells using a k-means clustering algorithm if `NUMBER OF INTEGRATION CELLS` is less than the number of sub-tet hex elements for an integration cell. Each sub-tet hex element will be used as a sub-cell if `NUMBER OF INTEGRATION CELLS` is equal to or greater than the available elements for an integration cell. The volumetric portion of the deformation gradient is calculated over the cell. These are the nodal or tetrahedral domains for the `NODAL_FROM_THEX` and `TET_FROM_THEX` options, respectively.

The `THEX ACTIVE NODES` command controls which nodes should be active and defaults to `NODES`. This command is only valid with `FORMULATION = NODAL_FROM_THEX` or `TET_FROM_THEX` where the input mesh has been subdivided using the `thex` command in cubit. The `thex` command splits every tetrahedral element into four hexahedral elements, requiring the creation of nodes at the edge and face midpoints and the centroid of the tetrahedron. Any of the nodes in the resulting mesh can be activated (have degrees-of-freedom associated with them). The option `NODES` activates the nodes that were in the original tet mesh. The options `EDGES`, `FACES` and `CENTROIDS` activates the nodes that the `thex` process created on the input tetrahedral edges, faces and centroids, respectively. This option is useful for specifying additional “bubble” DOFs when using advanced integration methods such as \(\mathbf{\overline{F}}\) [[2]].

## 14.2.3. Bond-based fracture

Warning

The bond-based fracture method is only compatible with the command `FORMULATION = NODAL_FROM_THEX`.

CRK currently supports an experimental fracture formulation known as bond-based fracture. This formulation creates a “bond” between each node-node pair and maintains a damage variable associated with the bond which degrades the incremental deformation gradient and resulting internal force [[3]]. The failure criterion is specified using the `FAILURE INDEX` command to specify the emph{0-based} state variable index of the material model in use along with `FAILURE CRITERION` command which specifies the threshold the state variable must exceed before the node is considered “damaged.” When both nodes of a bond have reached the criterion, the bond is marked as damaged and begins to degrade linearly over a distance given by the `BOND OPENING DISPLACEMENT` command. This bond opening displacement is given by the relative displacement between the two nodes (in an abuse of terminology, the “stretch” of the bond).

The criterion above requires both nodes of a bond to be flagged as “damaged” before the bond can begin degrading. An alternative criterion may be specified via the command `BOND SECONDARY CRITERION = <real>`\(u_\textrm{second}\). For any bond that fails due to the first criterion, this option creates a failure plane passing through that bond’s midpoint normal to the bond. Any bond which shares a node with the first failed bond and intersects the failure plane may also fail due to a second criterion given by

\[u_{ab} - \bar{u}_{ab} \ge u_\textrm{second},\]

where \(u_{ab} = \lVert\mathbf{u}_b - \mathbf{u}_a\rVert\) is the relative displacement between nodes \(a\) and \(b\), and \(\bar{u}_{ab}\) is the relative displacement at the time one of the bond’s nodes becomes damaged due to the criterion given by `FAILURE CRITERION` and `FAILURE INDEX` described above. The bond will then begin to degrade over the `BOND OPENING DISPLACEMENT` as before.

Warning

The bonds may be viewed using the `SHOW BONDS = ON` command. However, this currently only works on a single processor.

## 14.2.4. Output

Note, CRK shape functions in general do not possess the Kronecker delta property, i.e.

\[N_I(x_J) \neq \delta_{IJ}.\]

As such, nodal variables such as `displacement` and `velocity` are not true displacements or velocities but rather the emph{generalized coefficients} of those fields. In addition, other quantities of interest will also be generalized (e.g. mass, force, etc.). These generalized values must be converted to their physical counterparts by interpolating these coefficients with their respective shape functions. These quantities are stored in the variables `physical_displacement` and `physical_velocity`. As a workaround, one can provide an alternative name in the results output to switch the meaning of `displacement` and `physical_displacement` as follows:

    # Within Region scope...
    BEGIN RESULTS OUTPUT my_output
      ...
      NODAL VARIABLES = physical_displacement as displacement
      NODAL VARIABLES = displacement as gen_displacement
      ...
    END RESULTS OUTPUT my_output

which will output the displacement in a form amenable to Paraview or Ensight.

Warning

If an analysis contains both FEM and CRK formulations, then the `displacement` and `velocity` fields will be generalized for the CRK blocks and physical for the FEM blocks. Since `physical_displacement` does not exist in the FEM portions, the workaround mentioned above will output zero physical displacements for those blocks.

The output variables available for CRK are given below in Table 14.1 and Table 14.2 below. Bond-based fracture specific variables are given in Table 14.3.

Table 14.1 Nodal Variables for the CRK Formulation
Name | Type | Comments
--- | --- | ---
`physical_displacement` | `Vector_3D` | The true displacement predicted by CRK
`physical_velocity` | `Vector_3D` | The true velocity predicted by CRK
`displacement` | `Vector_3D` | The generalized displacement predicted by CRK, i.e. the displacement coefficients
`velocity` | `Vector_3D` | The generalized velocity predicted by CRK, i.e. the velocity coefficients
`nodal_status` | `Integer` | Indicates whether a node is active (carrying a DOF) or inactive
`max_edge_length` | `Real` | Maximum edge length connected to the node of the original, non `thex` mesh
`local_support_size` | `Real` | Local support size radius when `SHAPE FUNCTION TYPE` is one of `SPHERICAL` or `SPHERICAL_SURFACE_KRONECKER`
`num_neighbors` | `Integer` | Number of active neighbors which have shape function value at the node
`crk_contact_mass` | `Real` | Generalized mass scattered to inactive nodes for Dash purposes

Table 14.2 Element Variables for the CRK Formulation
Name | Type | Comments
--- | --- | ---
`cell_id` | `Integer` | The ID of the CRK strain smoothing cell
`sub_cell_id` | `Integer` | The ID of the CRK strain smoothing sub-cell
`cell_volume` | `Real` | The volume of the smoothing cell
`sub_cell_volume` | `Real` | The volume of the smoothing sub-cell

Table 14.3 Output Variables for Bond-Based FeFpmodelnameracture
Name | Type | Comments |
--- | --- | --- | ---
`node_damage` | Nodal | `Integer` | Flag indicating if the failure criterion is met at a given node’s material point
`damage_fraction` | Element | `Real` | Fraction of damaged bonds over total bonds associated with an integration cell
`bond_degradation` | Element | `Real` | Scalar between 0 and 1 to indicate level of damage in the bond where 1 indicates no damage and 0 indicates fully damaged (only available when `SHOW BONDS = ON`)

[1]

Jiun-Shyan Chen, Cheng-Tang Wu, Sangpil Yoon, and Yang You. A stabilized conforming nodal integration for galerkin mesh-free methods. International Journal for Numerical Methods in Engineering, 50(2):435–466, 2001.

[2] (1,2)

Georgios Moutsanidis, Jacob Koester, Michael Tupek, Yuri Bazilevs, and Jiun-Shyan Chen. Treatment of near-incompressibility in meshfree and immersed-particle methods. Computational Particle Mechanics, 2019.

[3]

Jacob Koester, Michael Tupek, and Scott Mitchell. An agile design-to-simulation workflow using a new conforming moving least squares method. Technical Report SAND2019-11851, Sandia National Laboratories, Albuquerque, NM, 2019.