Add coupled electrostatics test problem (idaholab#26105)

Co-authored-by: karthichockalingam <karthikeyan.chockalingam@stfc.ac.uk>

test/tests/kernels/vector_fe/coupled_electrostatics.i

@@ -0,0 +1,127 @@
+# Test for DivField and GradField kernels and VectorDivPenaltyDirichletBC bcs.
+# This test uses Raviart-Thomas elements to solve a model div-grad problem
+# in H(div). The problem is simply a div-grad formulation, u = -grad(p), and
+# div(u) = f, of the standard Poisson problem div(grad(p)) = -f.
+# Manufactured solution: p = cos(k*x)*sin(k*y)*cos(k*z).
+
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 6
+    ny = 6
+    nz = 6
+    xmax =  1
+    ymax =  1
+    zmax =  1
+    xmin = -1
+    ymin = -1
+    zmin = -1
+    elem_type = HEX27
+  []
+[]
+
+[Variables]
+  [u]
+    family = RAVIART_THOMAS
+    order = FIRST
+  []
+  [p]
+    family = MONOMIAL
+    order = CONSTANT
+  []
+  [lambda]
+    family = SCALAR
+    order = FIRST
+  []
+[]
+
+[Kernels]
+  [coefficient]
+    type = VectorFunctionReaction
+    variable = u
+    sign = negative
+  []
+  [gradient]
+    type = GradField
+    variable = u
+    coupled_scalar_variable = p
+  []
+  [divergence]
+    type = DivField
+    variable = p
+    coupled_vector_variable = u
+  []
+  [forcing]
+    type = BodyForce
+    variable = p
+    function = f
+  []
+  [mean_zero_p]
+    type = ScalarLagrangeMultiplier
+    variable = p
+    lambda = lambda
+  []
+[]
+
+[ScalarKernels]
+  [constraint]
+    type = AverageValueConstraint
+    variable = lambda
+    pp_name = pp
+    value = 0.0
+  []
+[]
+
+[BCs]
+  [sides]
+    type = VectorDivPenaltyDirichletBC
+    variable = u
+    function = s
+    penalty = 1e8
+    boundary = 'top bottom left right front back'
+  []
+[]
+
+[Functions]
+  [f]
+    type = ParsedFunction
+    expression = ${Mesh/gmg/dim}*k*k*cos(k*x)*sin(k*y)*cos(k*z)
+    symbol_names = k
+    symbol_values = 1.570796326794897 # pi/2
+  []
+  [s]
+    type = ParsedVectorFunction
+    expression_x =  k*sin(k*x)*sin(k*y)*cos(k*z)
+    expression_y = -k*cos(k*x)*cos(k*y)*cos(k*z)
+    expression_z =  k*cos(k*x)*sin(k*y)*sin(k*z)
+    symbol_names = k
+    symbol_values = 1.570796326794897 # pi/2
+  []
+[]
+
+[Postprocessors]
+  [pp]
+    type = ElementIntegralVariablePostprocessor
+    variable = p
+    execute_on = linear
+  []
+[]
+
+[Preconditioning]
+  [SMP]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = -pc_type
+  petsc_options_value = jacobi
+[]
+
+[Outputs]
+  exodus = true
+[]

test/tests/kernels/vector_fe/tests

@@ -151,4 +151,67 @@
     issues = "#13309"
     requirement = "The system shall error if the 'function' and 'function_x' parameters are both set within the ADVectorFunctionDirichletBC object."
   []
+  [coupled_electrostatics]
+    design = 'DivField.md GradField.md VectorDivPenaltyDirichletBC.md'
+    issues = '#26105'
+    [3d]
+      requirement = 'The system shall be capable of modeling a div-grad problem '
+                    'in H(div) using the following 3d Raviart-Thomas elements: '
+      [hex27]
+        type = Exodiff
+        input = coupled_electrostatics.i
+        exodiff = coupled_electrostatics_hex27.e
+        cli_args = 'Outputs/file_base=coupled_electrostatics_hex27'
+        detail = 'HEX27, and '
+      []
+      [tet14]
+        type = Exodiff
+        input = coupled_electrostatics.i
+        exodiff = coupled_electrostatics_tet14.e
+        cli_args = 'Mesh/gmg/elem_type=TET14 '
+                  'Outputs/file_base=coupled_electrostatics_tet14'
+        detail = 'TET14.'
+      []
+    []
+    [2d]
+      requirement = 'The system shall be capable of modeling a div-grad problem '
+                    'in H(div) using the following 2d Raviart-Thomas elements: '
+      [tri6]
+        type = Exodiff
+        input = coupled_electrostatics.i
+        exodiff = coupled_electrostatics_tri6.e
+        cli_args = 'Mesh/gmg/dim=2 Mesh/gmg/elem_type=TRI6 '
+                  'BCs/sides/boundary="top bottom left right" '
+                  'Outputs/file_base=coupled_electrostatics_tri6'
+        detail = 'TRI6, '
+      []
+      [tri7]
+        type = Exodiff
+        input = coupled_electrostatics.i
+        exodiff = coupled_electrostatics_tri7.e
+        cli_args = 'Mesh/gmg/dim=2 Mesh/gmg/elem_type=TRI7 '
+                  'BCs/sides/boundary="top bottom left right" '
+                  'Outputs/file_base=coupled_electrostatics_tri7'
+        detail = 'TRI7, '
+      []
+      [quad8]
+        type = Exodiff
+        input = coupled_electrostatics.i
+        exodiff = coupled_electrostatics_quad8.e
+        cli_args = 'Mesh/gmg/dim=2 Mesh/gmg/elem_type=QUAD8 '
+                  'BCs/sides/boundary="top bottom left right" '
+                  'Outputs/file_base=coupled_electrostatics_quad8'
+        detail = 'QUAD8, and '
+      []
+      [quad9]
+        type = Exodiff
+        input = coupled_electrostatics.i
+        exodiff = coupled_electrostatics_quad9.e
+        cli_args = 'Mesh/gmg/dim=2 Mesh/gmg/elem_type=QUAD9 '
+                  'BCs/sides/boundary="top bottom left right" '
+                  'Outputs/file_base=coupled_electrostatics_quad9'
+        detail = 'QUAD9.'
+      []
+    []
+  []
 []