No integration by parts on p

include/kernels/ADStokesStressDivergence.h

@@ -22,5 +22,5 @@ class ADStokesStressDivergence : public ADVectorKernel
   virtual ADReal computeQpResidual() override;
 
   const ADMaterialProperty<RankTwoTensor> & _stress;
-  const ADVariableValue & _pressure;
-};
+  const ADVariableGradient & _pressure;
+};

src/kernels/ADStokesStressDivergence.C

@@ -25,12 +25,12 @@ ADStokesStressDivergence::validParams()
 ADStokesStressDivergence::ADStokesStressDivergence(const InputParameters & parameters)
   : ADVectorKernel(parameters),
     _stress(getADMaterialPropertyByName<RankTwoTensor>("stress")),
-    _pressure(adCoupledValue("pressure"))
+    _pressure(adCoupledGradient("pressure"))
 {
 }
 
 ADReal
 ADStokesStressDivergence::computeQpResidual()
 {
-  return _stress[_qp].contract(_grad_test[_i][_qp]) - _grad_test[_i][_qp].tr() * _pressure[_qp];
+  return _stress[_qp].contract(_grad_test[_i][_qp]) + _test[_i][_qp] * _pressure[_qp];
 }

tests/stokes/manufactured.i

@@ -8,9 +8,9 @@
     ymax = 1
     zmin = 0
     zmax = 1
-    nx = 2
-    ny = 2
-    nz = 2
+    nx = 5
+    ny = 5
+    nz = 5
     elem_type = HEX20
   []
 []

tests/stokes/structural.i

@@ -0,0 +1,159 @@
+[Mesh]
+  second_order = true
+  [mesh]
+    type = CartesianMeshGenerator
+    dim = 3
+    dx = '1 1'
+    dy = '1 1'
+    dz = '2'
+    ix = '1 1'
+    iy = '1 1'
+    iz = '2'
+    subdomain_id = '1 2 3 4'
+  []
+[]
+
+[Variables]
+  [p]
+    order = FIRST
+    family = MONOMIAL
+  []
+  [u]
+    order = SECOND
+    family = LAGRANGE_VEC
+  []
+[]
+
+[AuxVariables]
+    [stress_xx]
+        family = MONOMIAL
+        order = CONSTANT
+    []
+    [stress_yy]
+        family = MONOMIAL
+        order = CONSTANT
+    []
+    [stress_zz]
+        family = MONOMIAL
+        order = CONSTANT
+    []
+[]
+
+[AuxKernels]
+    [stress_xx]
+        type = ADMaterialRankTwoTensorAux
+        variable = stress_xx
+        property = 'stress'
+        i = 0
+        j = 0
+    []
+    [stress_yy]
+        type = ADMaterialRankTwoTensorAux
+        variable = stress_yy
+        property = 'stress'
+        i = 1
+        j = 1
+    []
+    [stress_zz]
+        type = ADMaterialRankTwoTensorAux
+        variable = stress_zz
+        property = 'stress'
+        i = 2
+        j = 2
+    []
+[]
+
+[ICs]
+  [u]
+    type = VectorConstantIC
+    variable = u
+    x_value = 1e-15
+    y_value = 1e-15
+    z_value = 1e-15
+  []
+[]
+
+[Kernels]
+  [equil]
+    type = ADStokesStressDivergence
+    variable = u
+    pressure = p
+  []
+  [incompressible]
+    type = ADStokesIncompressibility
+    variable = p
+    velocity = u
+  []
+[]
+
+[Materials]
+  [strain]
+    type = StokesStrainRate
+    velocity = u
+  []
+  [stress]
+    type = StokesLinearViscous
+    mu = 1.0
+  []
+[]
+
+[BCs]
+  [fix_x]
+    type = ADVectorFunctionDirichletBC
+    variable = u
+    boundary = 'left'
+    set_y_comp = False
+    set_z_comp = False
+  []
+
+  [fix_y]
+    type = ADVectorFunctionDirichletBC
+    variable = u
+    boundary = 'bottom'
+    set_x_comp = False
+    set_z_comp = False
+  []
+
+  [fix_z]
+    type = ADVectorFunctionDirichletBC
+    variable = u
+    boundary = 'back'
+    set_x_comp = False
+    set_y_comp = False
+  []
+
+  [fix_p1]
+    type = ADDirichletBC
+    variable = p
+    boundary = 'front'
+    value = -1.0
+  []
+
+  [test]
+    type = ADVectorFunctionNeumannBC
+    boundary = 'front'
+    variable = u
+    function_x = '0'
+    function_y = '0'
+    function_z = '2.0/3.0'
+  []
+[]
+
+[Executioner]
+  type = Steady
+
+  solve_type = 'newton'
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'svd'
+  line_search = none
+
+  l_max_its = 10
+  l_tol = 1e-8
+  nl_max_its = 15
+  nl_rel_tol = 1e-10
+  nl_abs_tol = 1e-12
+[]
+
+[Outputs]
+  exodus = true
+[]