Merge pull request idaholab#15024 from lindsayad/ik-sidedness-for-gmsh

Allow interface kernel execution on gmsh sidesets

framework/src/interfacekernels/InterfaceKernel.C

@@ -197,6 +197,19 @@ template <typename T>
 void
 InterfaceKernelTempl<T>::computeResidual()
 {
+  // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
+  // associated only with the lower-dimensional geometric entity that is the boundary between two
+  // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
+  // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
+  // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
+  // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
+  // to execute in the former case. In the future we should remove this and add some kind of "block"
+  // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
+  // this return
+  if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
+      !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
+    return;
+
   // Compute the residual for this element
   computeElemNeighResidual(Moose::Element);
 
@@ -273,6 +286,19 @@ template <typename T>
 void
 InterfaceKernelTempl<T>::computeJacobian()
 {
+  // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
+  // associated only with the lower-dimensional geometric entity that is the boundary between two
+  // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
+  // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
+  // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
+  // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
+  // to execute in the former case. In the future we should remove this and add some kind of "block"
+  // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
+  // this return
+  if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
+      !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
+    return;
+
   computeElemNeighJacobian(Moose::ElementElement);
   computeElemNeighJacobian(Moose::NeighborNeighbor);
 }
@@ -313,6 +339,19 @@ template <typename T>
 void
 InterfaceKernelTempl<T>::computeElementOffDiagJacobian(unsigned int jvar)
 {
+  // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
+  // associated only with the lower-dimensional geometric entity that is the boundary between two
+  // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
+  // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
+  // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
+  // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
+  // to execute in the former case. In the future we should remove this and add some kind of "block"
+  // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
+  // this return
+  if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
+      !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
+    return;
+
   bool is_jvar_not_interface_var = true;
   if (jvar == _var.number())
   {
@@ -336,6 +375,19 @@ template <typename T>
 void
 InterfaceKernelTempl<T>::computeNeighborOffDiagJacobian(unsigned int jvar)
 {
+  // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
+  // associated only with the lower-dimensional geometric entity that is the boundary between two
+  // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
+  // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
+  // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
+  // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
+  // to execute in the former case. In the future we should remove this and add some kind of "block"
+  // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
+  // this return
+  if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
+      !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
+    return;
+
   bool is_jvar_not_interface_var = true;
   if (jvar == _var.number())
   {

test/tests/interfacekernels/gmsh_sidesets/coupled_value_coupled_flux.i

@@ -0,0 +1,86 @@
+[Mesh]
+  file = gmsh_mesh.msh
+[]
+
+[Variables]
+  [./u]
+    block = 6
+  [../]
+
+  [./v]
+    block = 5
+  [../]
+[]
+
+[Kernels]
+  [./diff_u]
+    type = CoeffParamDiffusion
+    variable = u
+    D = 4
+    block = 6
+  [../]
+  [./diff_v]
+    type = CoeffParamDiffusion
+    variable = v
+    D = 2
+    block = 5
+  [../]
+  [./source_u]
+    type = BodyForce
+    variable = u
+    value = 1
+  [../]
+[]
+
+[InterfaceKernels]
+ [./interface]
+   type = PenaltyInterfaceDiffusion
+   variable = u
+   neighbor_var = v
+   boundary = '1 2'
+   penalty = 1e6
+ [../]
+[]
+
+[BCs]
+  [./u]
+    type = VacuumBC
+    variable = u
+    boundary = 4
+  [../]
+  [./v]
+    type = VacuumBC
+    variable = v
+    boundary = 3
+  [../]
+[]
+
+[Postprocessors]
+  [./u_int]
+    type = ElementIntegralVariablePostprocessor
+    variable = u
+    block = 6
+  [../]
+  [./v_int]
+    type = ElementIntegralVariablePostprocessor
+    variable = v
+    block = 5
+  [../]
+[]
+
+[Preconditioning]
+  [./smp]
+    type = SMP
+    full = true
+  [../]
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+[]
+
+[Outputs]
+  exodus = true
+  print_linear_residuals = true
+[]

test/tests/interfacekernels/gmsh_sidesets/gmsh_mesh.msh

@@ -0,0 +1,57 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$PhysicalNames
+6
+1 1 "Interface_top"
+1 2 "Interface_bottom"
+1 3 "Boundary_bottom_left"
+1 4 "Boundary_top_right"
+2 5 "block_1"
+2 6 "block_2"
+$EndPhysicalNames
+$Nodes
+13
+1 0 0 0
+2 0 2 0
+3 2 2 0
+4 2 0 0
+5 1 0 0
+6 1 1 0
+7 0 1 0
+8 2 0.999999999997388 0
+9 1.000000000004118 2 0
+10 0.5 0.5 0
+11 1.50000000000103 1.499999999999347 0
+12 0.5000000000010296 1.5 0
+13 1.5 0.499999999998694 0
+$EndNodes
+$Elements
+26
+1 1 2 3 1 1 5
+2 1 2 4 2 5 4
+3 1 2 4 3 4 8
+4 1 2 4 3 8 3
+5 1 2 4 4 3 9
+6 1 2 4 4 9 2
+7 1 2 4 5 2 7
+8 1 2 3 6 7 1
+9 1 2 2 7 5 6
+10 1 2 1 8 6 7
+11 2 2 5 1 1 5 10
+12 2 2 5 1 7 1 10
+13 2 2 5 1 5 6 10
+14 2 2 5 1 6 7 10
+15 2 2 6 2 8 3 11
+16 2 2 6 2 9 2 12
+17 2 2 6 2 9 6 11
+18 2 2 6 2 6 8 11
+19 2 2 6 2 6 9 12
+20 2 2 6 2 5 4 13
+21 2 2 6 2 8 6 13
+22 2 2 6 2 6 5 13
+23 2 2 6 2 2 7 12
+24 2 2 6 2 7 6 12
+25 2 2 6 2 3 9 11
+26 2 2 6 2 4 8 13
+$EndElements

test/tests/interfacekernels/gmsh_sidesets/tests

@@ -0,0 +1,10 @@
+[Tests]
+  [./test]
+    type = 'Exodiff'
+    input = 'coupled_value_coupled_flux.i'
+    exodiff = 'coupled_value_coupled_flux_out.e'
+    issues = "#15203"
+    design = "InterfaceKernels/index.md"
+    requirement = "The system shall be able to execute interfacial conditions on sidesets generated by gmsh."
+  [../]
+[]