diff --git a/modules/combined/test/tests/phase_field_fracture/crack2d_no_split.i b/modules/combined/test/tests/phase_field_fracture/crack2d_no_split.i
new file mode 100644
index 000000000000..6a5fcd11039a
--- /dev/null
+++ b/modules/combined/test/tests/phase_field_fracture/crack2d_no_split.i
@@ -0,0 +1,183 @@
+#This input uses PhaseField-Nonconserved Action to add phase field fracture bulk rate kernels
+[Mesh]
+  [gen]
+    type = GeneratedMeshGenerator
+    dim = 2
+    nx = 20
+    ny = 10
+    ymax = 0.5
+  []
+  [./noncrack]
+    type = BoundingBoxNodeSetGenerator
+    new_boundary = noncrack
+    bottom_left = '0.5 0 0'
+    top_right = '1 0 0'
+    input = gen
+  [../]
+[]
+
+[GlobalParams]
+  displacements = 'disp_x disp_y'
+[]
+
+[Modules]
+  [./PhaseField]
+    [./Nonconserved]
+      [./c]
+        free_energy = F
+        kappa = kappa_op
+        mobility = L
+      [../]
+    [../]
+  [../]
+  [./TensorMechanics]
+    [./Master]
+      [./mech]
+        add_variables = true
+        strain = SMALL
+        additional_generate_output = 'stress_yy'
+        save_in = 'resid_x resid_y'
+      [../]
+    [../]
+  [../]
+[]
+
+[AuxVariables]
+  [./resid_x]
+  [../]
+  [./resid_y]
+  [../]
+[]
+
+[Kernels]
+  [./solid_x]
+    type = PhaseFieldFractureMechanicsOffDiag
+    variable = disp_x
+    component = 0
+    c = c
+  [../]
+  [./solid_y]
+    type = PhaseFieldFractureMechanicsOffDiag
+    variable = disp_y
+    component = 1
+    c = c
+  [../]
+[]
+
+[BCs]
+  [./ydisp]
+    type = FunctionDirichletBC
+    variable = disp_y
+    boundary = top
+    function = 't'
+  [../]
+  [./yfix]
+    type = DirichletBC
+    variable = disp_y
+    boundary = noncrack
+    value = 0
+  [../]
+  [./xfix]
+    type = DirichletBC
+    variable = disp_x
+    boundary = top
+    value = 0
+  [../]
+[]
+
+[Materials]
+  [./pfbulkmat]
+    type = GenericConstantMaterial
+    prop_names = 'gc_prop l visco'
+    prop_values = '1e-3 0.04 1e-4'
+  [../]
+  [./define_mobility]
+    type = ParsedMaterial
+    material_property_names = 'gc_prop visco'
+    f_name = L
+    function = '1.0/(gc_prop * visco)'
+  [../]
+  [./define_kappa]
+    type = ParsedMaterial
+    material_property_names = 'gc_prop l'
+    f_name = kappa_op
+    function = 'gc_prop * l'
+  [../]
+  [./elasticity_tensor]
+    type = ComputeElasticityTensor
+    C_ijkl = '120.0 80.0'
+    fill_method = symmetric_isotropic
+  [../]
+  [./damage_stress]
+    type = ComputeLinearElasticPFFractureStress
+    c = c
+    E_name = 'elastic_energy'
+    D_name = 'degradation'
+    F_name = 'local_fracture_energy'
+    decomposition_type = none
+  [../]
+  [./degradation]
+    type = DerivativeParsedMaterial
+    f_name = degradation
+    args = 'c'
+    function = '(1.0-c)^2*(1.0 - eta) + eta'
+    constant_names       = 'eta'
+    constant_expressions = '0.0'
+    derivative_order = 2
+  [../]
+  [./local_fracture_energy]
+    type = DerivativeParsedMaterial
+    f_name = local_fracture_energy
+    args = 'c'
+    material_property_names = 'gc_prop l'
+    function = 'c^2 * gc_prop / 2 / l'
+    derivative_order = 2
+  [../]
+  [./fracture_driving_energy]
+    type = DerivativeSumMaterial
+    args = c
+    sum_materials = 'elastic_energy local_fracture_energy'
+    derivative_order = 2
+    f_name = F
+  [../]
+[]
+
+[Postprocessors]
+  [./resid_x]
+    type = NodalSum
+    variable = resid_x
+    boundary = 2
+  [../]
+  [./resid_y]
+    type = NodalSum
+    variable = resid_y
+    boundary = 2
+  [../]
+[]
+
+[Preconditioning]
+  [./smp]
+    type = SMP
+    full = true
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+
+  solve_type = PJFNK
+  petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_ksp_type -sub_pc_type -pc_asm_overlap'
+  petsc_options_value = 'asm      31                  preonly       lu           1'
+
+  nl_rel_tol = 1e-8
+  l_max_its = 10
+  nl_max_its = 10
+
+  dt = 1e-4
+  dtmin = 1e-4
+  num_steps = 2
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/modules/combined/test/tests/phase_field_fracture/gold/crack2d_no_split_out.e b/modules/combined/test/tests/phase_field_fracture/gold/crack2d_no_split_out.e
new file mode 100644
index 000000000000..cd2f7c1a52ca
Binary files /dev/null and b/modules/combined/test/tests/phase_field_fracture/gold/crack2d_no_split_out.e differ
diff --git a/modules/combined/test/tests/phase_field_fracture/tests b/modules/combined/test/tests/phase_field_fracture/tests
index aeb5e848280c..319ddf1b2e5f 100644
--- a/modules/combined/test/tests/phase_field_fracture/tests
+++ b/modules/combined/test/tests/phase_field_fracture/tests
@@ -19,6 +19,16 @@
     issues = '#12868'
     requirement = 'The phase-field fracture model shall predict damage evolution using strain volumetric/deviatoric decomposition approach for an isotropic elastic material.'
   [../]
+  [./crack2d_no_split]
+    type = 'Exodiff'
+    input = 'crack2d_no_split.i'
+    exodiff = 'crack2d_no_split_out.e'
+    abs_zero = 1e-09
+    use_old_floor = True
+    design = 'ComputeLinearElasticPFFractureStress.md'
+    issues = '#12868 #15695'
+    requirement = 'The phase-field fracture model shall predict damage evolution with no decomposition approach for an isotropic elastic material.'
+  [../]
   [./crack2d_linear_fracture_energy]
     type = 'Exodiff'
     input = 'crack2d_linear_fracture_energy.i'
diff --git a/modules/tensor_mechanics/src/materials/ComputeLinearElasticPFFractureStress.C b/modules/tensor_mechanics/src/materials/ComputeLinearElasticPFFractureStress.C
index e874be9e491b..6b77177c42c6 100644
--- a/modules/tensor_mechanics/src/materials/ComputeLinearElasticPFFractureStress.C
+++ b/modules/tensor_mechanics/src/materials/ComputeLinearElasticPFFractureStress.C
@@ -211,14 +211,14 @@ ComputeLinearElasticPFFractureStress::computeQpStress()
       break;
     default:
     {
-      _stress[_qp] = _D[_qp] * _elasticity_tensor[_qp] * _mechanical_strain[_qp] -
-                     _pressure[_qp] * I2 * _I[_qp];
-      F_pos = (_stress[_qp]).doubleContraction(_mechanical_strain[_qp]) / 2.0;
+      RankTwoTensor stress = _elasticity_tensor[_qp] * _mechanical_strain[_qp];
+      F_pos = stress.doubleContraction(_mechanical_strain[_qp]) / 2.0;
       F_neg = 0.0;
       if (_use_current_hist)
-        _d2Fdcdstrain[_qp] = _stress[_qp] * _dDdc[_qp];
+        _d2Fdcdstrain[_qp] = stress * _dDdc[_qp];
 
-      _dstress_dc[_qp] = _stress[_qp] * _dDdc[_qp] - _pressure[_qp] * I2 * _dIdc[_qp];
+      _stress[_qp] = _D[_qp] * stress - _pressure[_qp] * I2 * _I[_qp];
+      _dstress_dc[_qp] = stress * _dDdc[_qp] - _pressure[_qp] * I2 * _dIdc[_qp];
       _Jacobian_mult[_qp] = _D[_qp] * _elasticity_tensor[_qp];
     }
   }