Merge pull request #25838 from farscape-project/tet14

Add support for TET14 geometric elements

framework/include/mesh/MooseMesh.h

@@ -59,7 +59,7 @@ const std::string LIST_GEOM_ELEM = "EDGE EDGE2 EDGE3 EDGE4 "
                                    "QUAD QUAD4 QUAD8 QUAD9 "
                                    "TRI TRI3 TRI6 TRI7 "
                                    "HEX HEX8 HEX20 HEX27 "
-                                   "TET TET4 TET10 "
+                                   "TET TET4 TET10 TET14 "
                                    "PRISM PRISM6 PRISM15 PRISM18 "
                                    "PYRAMID PYRAMID5 PYRAMID13 PYRAMID14";
 

modules/heat_transfer/test/tests/directional_flux_bc/gold/3d_tet14.csv

@@ -0,0 +1,3 @@
+time,light
+0,0
+0.01,0.16229642731081

modules/heat_transfer/test/tests/directional_flux_bc/tests

@@ -79,6 +79,14 @@
       detail = 'TET10'
       heavy = true
     []
+    [tet14]
+      type = CSVDiff
+      input = 3d_elem.i
+      csvdiff = 3d_tet14.csv
+      cli_args = 'GlobalParams/elem_type=tet14 GlobalParams/order=third Outputs/file_base=3d_tet14'
+      detail = 'TET14'
+      heavy = true
+    []
 
     [hex8]
       type = CSVDiff

modules/ray_tracing/src/raytracing/TraceRay.C

@@ -19,6 +19,7 @@
 // libMesh includes
 #include "libmesh/cell_tet4.h"
 #include "libmesh/cell_tet10.h"
+#include "libmesh/cell_tet14.h"
 #include "libmesh/cell_hex8.h"
 #include "libmesh/cell_hex20.h"
 #include "libmesh/cell_hex27.h"
@@ -220,6 +221,15 @@ TraceRay::exitsElem(const Elem * elem,
                                            intersection_distance,
                                            normals);
       break;
+    case TET14:
+      intersected = exitsElem<Tet14, Tet4>(elem,
+                                           incoming_side,
+                                           intersection_point,
+                                           intersected_side,
+                                           intersected_extrema,
+                                           intersection_distance,
+                                           normals);
+      break;
     case PYRAMID13:
       intersected = exitsElem<Pyramid13, Pyramid5>(elem,
                                                    incoming_side,

modules/ray_tracing/src/raytracing/TraceRayTools.C

@@ -25,6 +25,7 @@
 #include "libmesh/cell_pyramid14.h"
 #include "libmesh/cell_tet4.h"
 #include "libmesh/cell_tet10.h"
+#include "libmesh/cell_tet14.h"
 #include "libmesh/edge_edge2.h"
 #include "libmesh/enum_to_string.h"
 #include "libmesh/face_quad4.h"
@@ -37,7 +38,7 @@ namespace TraceRayTools
 {
 
 const std::set<int> TRACEABLE_ELEMTYPES = {
-    HEX8, HEX20, HEX27, QUAD4, QUAD8,    QUAD9,     TET4,      TET10,  TRI3,    TRI6,
+    HEX8, HEX20, HEX27, QUAD4, QUAD8,    QUAD9,     TET4,      TET10,  TET14,   TRI3,   TRI6,
     TRI7, EDGE2, EDGE3, EDGE4, PYRAMID5, PYRAMID13, PYRAMID14, PRISM6, PRISM15, PRISM18};
 const std::set<int> ADAPTIVITY_TRACEABLE_ELEMTYPES = {QUAD4, HEX8, TRI3, TET4, EDGE2};
 
@@ -291,6 +292,9 @@ findEdgeNeighbors(
       case TET10:
         return findEdgeNeighborsWithinEdgeInternal<Tet10>(
             candidate, elem, node1, node2, edge_length, info);
+      case TET14:
+        return findEdgeNeighborsWithinEdgeInternal<Tet14>(
+            candidate, elem, node1, node2, edge_length, info);
       case PYRAMID13:
         return findEdgeNeighborsWithinEdgeInternal<Pyramid13>(
             candidate, elem, node1, node2, edge_length, info);
@@ -597,6 +601,7 @@ withinEdge(const Elem * elem,
       return withinEdgeTempl<Hex8>(elem, point, extrema, tolerance);
     case TET4:
     case TET10:
+    case TET14:
       return withinEdgeTempl<Tet4>(elem, point, extrema, tolerance);
     case PYRAMID5:
     case PYRAMID13:
@@ -634,6 +639,7 @@ atVertexOnSide(const Elem * elem, const Point & point, const unsigned short side
       return atVertexOnSideTempl<Tri3>(elem, point, side);
     case TET4:
     case TET10:
+    case TET14:
       return atVertexOnSideTempl<Tet4>(elem, point, side);
     case PYRAMID5:
     case PYRAMID13:
@@ -699,6 +705,7 @@ withinEdgeOnSide(const Elem * const elem,
       return withinEdgeOnSideTempl<Hex8>(elem, point, side, extrema);
     case TET4:
     case TET10:
+    case TET14:
       return withinEdgeOnSideTempl<Tet4>(elem, point, side, extrema);
     case PYRAMID5:
     case PYRAMID13:

modules/ray_tracing/test/tests/traceray/tet14/gold/centroid_aq_out.csv

@@ -0,0 +1,3 @@
+time,distance_difference,expected_distance,total_distance
+0,0,0,0
+1,0,4634.7221896827,4634.7221896827

modules/ray_tracing/test/tests/traceray/tet14/gold/ctc_out.csv

@@ -0,0 +1,3 @@
+time,distance_difference,expected_distance,total_distance
+0,0,0,0
+1,0,193350.92635802,193350.92635802

modules/ray_tracing/test/tests/traceray/tet14/gold/ctv_out.csv

@@ -0,0 +1,3 @@
+time,distance_difference,expected_distance,total_distance
+0,0,0,0
+1,3.4106051316485e-13,1009.5543571299,1009.5543571299

modules/ray_tracing/test/tests/traceray/tet14/gold/ete_out.csv

@@ -0,0 +1,3 @@
+time,distance_difference,expected_distance,total_distance
+0,0,0,0
+1,0,11020.49719605,11020.49719605

modules/ray_tracing/test/tests/traceray/tet14/gold/side_aq_out.csv

@@ -0,0 +1,3 @@
+time,distance_difference,expected_distance,total_distance
+0,0,0,0
+1,0,7002.2125376021,7002.2125376021

modules/ray_tracing/test/tests/traceray/tet14/gold/vtv_out.csv

@@ -0,0 +1,3 @@
+time,distance_difference,expected_distance,total_distance
+0,0,0,0
+1,0,3019.8969044772,3019.8969044772

modules/ray_tracing/test/tests/traceray/tet14/tests

@@ -0,0 +1,121 @@
+[Tests]
+  design = 'RayTracingStudy.md'
+  issues = '#16028'
+
+  [tests]
+    requirement = 'The system shall be able to trace rays within planar TET14 elements '
+    [vtv]
+      type = 'CSVDiff'
+      input = '../lots.i'
+      csvdiff = 'vtv_out.csv'
+      cli_args = 'Mesh/gmg/dim=3
+                  Mesh/gmg/nx=3
+                  Mesh/gmg/ny=3
+                  Mesh/gmg/nz=3
+                  Mesh/gmg/elem_type=tet14
+                  UserObjects/lots/vertex_to_vertex=true
+                  Outputs/file_base=vtv_out
+                  RayBCs/active=kill_3d'
+      allow_test_objects = true
+      abs_zero = 1e-8
+      method = '!dbg'
+
+      detail = 'from boundary vertices on boundary sides in the direction of the vertices on the other side of the same boundary element, '
+    []
+
+    [ctv]
+      type = 'CSVDiff'
+      input = '../lots.i'
+      csvdiff = 'ctv_out.csv'
+      cli_args = 'Mesh/gmg/dim=3
+                  Mesh/gmg/nx=3
+                  Mesh/gmg/ny=3
+                  Mesh/gmg/nz=3
+                  Mesh/gmg/elem_type=tet14
+                  UserObjects/lots/centroid_to_vertex=true
+                  Outputs/file_base=ctv_out
+                  RayBCs/active=kill_3d'
+      allow_test_objects = true
+      abs_zero = 1e-8
+      method = '!dbg'
+
+      detail = 'from the centroids of boundary sides in the direction of all other vertices in the same boundary element, '
+    []
+
+    [ctc]
+      type = 'CSVDiff'
+      input = '../lots.i'
+      csvdiff = 'ctc_out.csv'
+      cli_args = 'Mesh/gmg/dim=3
+                  Mesh/gmg/nx=3
+                  Mesh/gmg/ny=3
+                  Mesh/gmg/nz=3
+                  Mesh/gmg/elem_type=tet14
+                  UserObjects/lots/centroid_to_centroid=true
+                  Outputs/file_base=ctc_out
+                  RayBCs/active=kill_3d'
+      allow_test_objects = true
+      abs_zero = 1e-8
+      method = '!dbg'
+
+      detail = 'from the centroids of all boundary sides to the centroids of all other boundary elements, '
+    []
+
+    [ete]
+      type = 'CSVDiff'
+      input = '../lots.i'
+      csvdiff = 'ete_out.csv'
+      cli_args = 'Mesh/gmg/dim=3
+                  Mesh/gmg/nx=3
+                  Mesh/gmg/ny=3
+                  Mesh/gmg/nz=3
+                  Mesh/gmg/elem_type=tet14
+                  UserObjects/lots/edge_to_edge=true
+                  Outputs/file_base=ete_out
+                  RayBCs/active=kill_3d'
+      allow_test_objects = true
+      abs_zero = 1e-8
+      method = '!dbg'
+
+      detail = 'from the centroids of boundary edges in the direction of all other edge centroids in the same element, '
+    []
+
+    [side_aq]
+      type = 'CSVDiff'
+      input = '../lots.i'
+      csvdiff = 'side_aq_out.csv'
+      cli_args = 'Mesh/gmg/dim=3
+                  Mesh/gmg/nx=3
+                  Mesh/gmg/ny=3
+                  Mesh/gmg/nz=3
+                  Mesh/gmg/elem_type=tet14
+                  UserObjects/lots/side_aq=true
+                  Outputs/file_base=side_aq_out
+                  RayBCs/active=kill_3d'
+      allow_test_objects = true
+      abs_zero = 1e-8
+      method = '!dbg'
+
+      detail = 'from boundary sides centroids using an angular quadrature, '
+    []
+
+    [centroid_aq]
+      type = 'CSVDiff'
+      input = '../lots.i'
+      csvdiff = 'centroid_aq_out.csv'
+      cli_args = 'Mesh/gmg/dim=3
+                  Mesh/gmg/nx=3
+                  Mesh/gmg/ny=3
+                  Mesh/gmg/nz=3
+                  Mesh/gmg/elem_type=tet14
+                  UserObjects/lots/centroid_aq=true
+                  Outputs/file_base=centroid_aq_out
+                  RayBCs/active=kill_3d'
+      allow_test_objects = true
+      abs_zero = 1e-8
+      method = '!dbg'
+
+      detail = 'and from boundary element centroids using an angular quadrature.'
+    []
+  []
+[]

modules/ray_tracing/unit/src/TraceRayToolsTest.C

@@ -43,6 +43,7 @@ traceRayToolsTestMesh(const int type)
     case HEX27:
     case TET4:
     case TET10:
+    case TET14:
     case PYRAMID5:
     case PYRAMID13:
     case PYRAMID14:

modules/xfem/src/efa/EFAElement3D.C

@@ -33,12 +33,14 @@ EFAElement3D::EFAElement3D(unsigned int eid, unsigned int n_nodes, unsigned int
   if (_num_faces == 4)
   {
     _num_vertices = 4;
-    if (_num_nodes == 10)
+    if (_num_nodes == 14)
+      _num_interior_face_nodes = 4;
+    else if (_num_nodes == 10)
       _num_interior_face_nodes = 3;
     else if (_num_nodes == 4)
       _num_interior_face_nodes = 0;
     else
-      EFAError("In EFAelement3D the supported TET element types are TET4 and TET10");
+      EFAError("In EFAelement3D the supported TET element types are TET4, TET10 and TET14");
   }
   else if (_num_faces == 6)
   {
@@ -53,7 +55,8 @@ EFAElement3D::EFAElement3D(unsigned int eid, unsigned int n_nodes, unsigned int
       EFAError("In EFAelement3D the supported HEX element types are HEX8, HEX20 and HEX27");
   }
   else
-    EFAError("In EFAelement3D the supported element types are TET4, TET10, HEX8, HEX20 and HEX27");
+    EFAError("In EFAelement3D the supported element types are TET4, TET10, TET14, HEX8, HEX20 and "
+             "HEX27");
   setLocalCoordinates();
 }
 
@@ -230,6 +233,13 @@ EFAElement3D::setLocalCoordinates()
                   1
 
     */
+    /*
+      TET14 elements also include four face nodes:
+      Node 10, centroid on side 0, arithmetic mean of 0/1/2 or 4/5/6
+      Node 11, centroid on side 1, arithmetic mean of 0/1/3 or 4/7/8
+      Node 12, centroid on side 2, arithmetic mean of 1/2/3 or 5/8/9
+      Node 13, centroid on side 3, arithmetic mean of 0/2/3 or 6/7/9
+    */
     _local_node_coor.resize(_num_nodes);
     _local_node_coor[0] = EFAPoint(0.0, 0.0, 0.0);
     _local_node_coor[1] = EFAPoint(1.0, 0.0, 0.0);
@@ -245,6 +255,14 @@ EFAElement3D::setLocalCoordinates()
       _local_node_coor[8] = EFAPoint(0.5, 0.0, 0.5);
       _local_node_coor[9] = EFAPoint(0.0, 0.5, 0.5);
     }
+
+    if (_num_nodes > 10)
+    {
+      _local_node_coor[10] = EFAPoint(1 / 3., 1 / 3., 0.0);
+      _local_node_coor[11] = EFAPoint(1 / 3., 0.0, 1 / 3.);
+      _local_node_coor[12] = EFAPoint(1 / 3., 1 / 3., 1 / 3.);
+      _local_node_coor[13] = EFAPoint(0.0, 1 / 3., 1 / 3.);
+    }
   }
   else
     EFAError("EFAElement3D: number of faces should be either 4(TET) or 6(HEX).");
@@ -1433,7 +1451,8 @@ EFAElement3D::createFaces()
                                               {10, 14, 18, 15, 23},
                                               {11, 15, 19, 12, 24},
                                               {16, 17, 18, 19, 25}};
-  int tet_interior_face_node_indices[4][3] = {{4, 5, 6}, {4, 7, 8}, {5, 8, 9}, {6, 7, 9}};
+  int tet_interior_face_node_indices[4][4] = {
+      {4, 5, 6, 10}, {4, 7, 8, 11}, {5, 8, 9, 12}, {6, 7, 9, 13}};
 
   _faces = std::vector<EFAFace *>(_num_faces, nullptr);
   if (_num_nodes == 8 || _num_nodes == 20 || _num_nodes == 27)
@@ -1450,7 +1469,7 @@ EFAElement3D::createFaces()
         _faces[i]->setInteriorFaceNode(k, _nodes[hex_interior_face_node_indices[i][k]]);
     }
   }
-  else if (_num_nodes == 4 || _num_nodes == 10)
+  else if (_num_nodes == 4 || _num_nodes == 10 || _num_nodes == 14)
   {
     if (_num_faces != 4)
       EFAError("num_faces of tets must be 4");

modules/xfem/src/efa/ElementFragmentAlgorithm.C

@@ -134,15 +134,9 @@ ElementFragmentAlgorithm::add3DElement(const std::vector<unsigned int> & quad, u
 {
   unsigned int num_nodes = quad.size();
   unsigned int num_faces = 0;
-  if (num_nodes == 27)
+  if (num_nodes == 8 || num_nodes == 20 || num_nodes == 27)
     num_faces = 6;
-  else if (num_nodes == 20)
-    num_faces = 6;
-  else if (num_nodes == 8)
-    num_faces = 6;
-  else if (num_nodes == 4)
-    num_faces = 4;
-  else if (num_nodes == 10)
+  else if (num_nodes == 4 || num_nodes == 10 || num_nodes == 14)
     num_faces = 4;
   else
     EFAError("In add3DElement element with id: ", id, " has invalid num_nodes");

modules/xfem/test/tests/high_order_elements/tests

@@ -54,6 +54,15 @@
     unique_id = true
     requirement = 'The XFEM module shall permit modeling of discontinuities represented with XFEM with second-order elements in 3D using TET10 elements when the XFEM cutting plane is prescribed using a level set field.'
   [../]
+  [./diffusion_3d_tet14]
+    type = Exodiff
+    input = diffusion_3d.i
+    exodiff = 'diffusion_3d_tet14_out.e'
+    cli_args = 'GlobalParams/order=THIRD Mesh/elem_type=TET14 Outputs/file_base=diffusion_3d_tet14_out'
+    map = false
+    unique_id = true
+    requirement = 'The XFEM module shall permit modeling of discontinuities represented with XFEM with third-order elements in 3D using TET14 elements when the XFEM cutting plane is prescribed using a level set field.'
+  [../]
   [./diffusion_3d_hex20]
     type = Exodiff
     input = diffusion_3d.i

test/tests/mesh/high_order_elems/tests

@@ -104,6 +104,15 @@
     requirement = 'The system shall support refinement of TET10 mesh elements.'
   [../]
 
+  [./test_tet14_refine]
+    type = 'Exodiff'
+    input = 'high_order_elems.i'
+    exodiff = 'high_order_elems_tet14_refine_out.e'
+    cli_args = 'Mesh/uniform_refine=1 Mesh/elem_type=TET14 Variables/u/order=THIRD Mesh/dim=3 Mesh/nx=1 Mesh/ny=1 Mesh/nz=1 Outputs/file_base=high_order_elems_tet14_refine_out'
+
+    requirement = 'The system shall support refinement of TET14 mesh elements.'
+  [../]
+
   [./test_prism6_refine]
     type = 'Exodiff'
     input = 'high_order_elems.i'

test/tests/meshgenerators/mesh_diagnostics_generator/tests

@@ -201,6 +201,15 @@
 
       detail = 'tetrahedral elements of second order,'
     []
+    [tet_third_order]
+      type = 'RunException'
+      input = 'detect_amr_tet.i'
+      mesh_mode = replicated
+      cli_args = 'Mesh/cut_one/elem_type=TET14 Mesh/diag/search_for_adaptivity_nonconformality=ERROR --mesh-only'
+      expect_err = 'Number of non-conformal nodes likely due to mesh refinement detected by heuristic: 8'
+
+      detail = 'tetrahedral elements of third order,'
+    []
   []
 
   [amr_detection_avoid_false_positives]