Merge topic 'tpl-verdict-1.4.2'

82b7c7a9 verdict: update spdx location and version to 1.4.2
0a745a92 Merge branch 'upstream-verdict' into tpl-verdict-1.4.2
e73a7f89 verdict 2024-12-30 (99f88172)
a410a262 verdict: update version to 1.4.2

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: buildbot <buildbot@kitware.com>
Reviewed-by: Cory Quammen <cory.quammen@kitware.com>
Merge-request: !11792

ThirdParty/verdict/CMakeLists.txt

@@ -7,9 +7,9 @@ vtk_module_third_party(
     SPDX_COPYRIGHT_TEXT
       "Copyright 2003,2006,2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS)."
     SPDX_DOWNLOAD_LOCATION
-      "git+https://gitlab.kitware.com/third-party/verdict.git@for/vtk-20220303-1.4.0"
+      "git+https://gitlab.kitware.com/third-party/verdict.git@for/vtk-20241230-1.4.2"
     VERSION
-      "1.4.0"
+      "1.4.2"
     STANDARD_INCLUDE_DIRS
   EXTERNAL
     PACKAGE Verdict

ThirdParty/verdict/update.sh

@@ -8,7 +8,7 @@ readonly name="verdict"
 readonly ownership="Verdict Upstream <kwrobot@kitware.com>"
 readonly subtree="ThirdParty/$name/vtk$name"
 readonly repo="https://gitlab.kitware.com/third-party/verdict.git"
-readonly tag="for/vtk-20231029-1.4.0"
+readonly tag="for/vtk-20241230-1.4.2"
 readonly paths="
 .gitattributes
 V_*

ThirdParty/verdict/vtkverdict/V_EdgeMetric.cpp

@@ -27,12 +27,30 @@ namespace VERDICT_NAMESPACE
 /*!\brief Length of and edge.
  * Length is calculated by taking the distance between the end nodes.
  */
-double edge_length(int /*num_nodes*/, const double coordinates[][3])
+double edge_length(int num_nodes, const double coordinates[][3])
 {
-
-  double x = coordinates[1][0] - coordinates[0][0];
-  double y = coordinates[1][1] - coordinates[0][1];
-  double z = coordinates[1][2] - coordinates[0][2];
-  return (double)(std::sqrt(x * x + y * y + z * z));
+  double edge_length = 0.0;
+
+  if (2 == num_nodes)
+  {  
+    double x = coordinates[1][0] - coordinates[0][0];
+    double y = coordinates[1][1] - coordinates[0][1];
+    double z = coordinates[1][2] - coordinates[0][2];
+    edge_length = (double)(std::sqrt(x * x + y * y + z * z));   
+  }
+  if (3 == num_nodes)
+  {
+    double x = coordinates[2][0] - coordinates[0][0];
+    double y = coordinates[2][1] - coordinates[0][1];
+    double z = coordinates[2][2] - coordinates[0][2];
+    edge_length += (double)(std::sqrt(x * x + y * y + z * z));
+
+    x = coordinates[2][0] - coordinates[1][0];
+    y = coordinates[2][1] - coordinates[1][1];
+    z = coordinates[2][2] - coordinates[1][2];
+    edge_length += (double)(std::sqrt(x * x + y * y + z * z));
+  }
+
+  return edge_length;
 }
 } // namespace verdict

ThirdParty/verdict/vtkverdict/V_HexMetric.cpp

@@ -105,7 +105,7 @@ static double compute_tet_volume(VerdictVector& v1, VerdictVector& v2, VerdictVe
 // Compute interior node
 VerdictVector hex20_auxillary_node_coordinate(const double coordinates[][3])
 {
-  VerdictVector aux_node(0, 0, 0);
+  VerdictVector aux_node(0.0, 0.0, 0.0);
   for (int i = 0; i < 8; i++)
   {
     VerdictVector tmp_vec(coordinates[i][0], coordinates[i][1], coordinates[i][2]);

ThirdParty/verdict/vtkverdict/V_QuadMetric.cpp

@@ -24,6 +24,7 @@
 #include "verdict_defines.hpp"
 
 #include <algorithm>
+#include <vector>
 
 namespace VERDICT_NAMESPACE
 {
@@ -459,16 +460,10 @@ double quad_aspect_ratio(int /*num_nodes*/, const double coordinates[][3])
   double mb = c1 > d1 ? c1 : d1;
   double hm = ma > mb ? ma : mb;
 
-  VerdictVector ab = edges[0] * edges[1];
-  VerdictVector cd = edges[2] * edges[3];
-  double denominator = ab.length() + cd.length();
-
-  if (denominator < VERDICT_DBL_MIN)
-  {
-    return (double)VERDICT_DBL_MAX;
-  }
+  double corner_areas[4];
+  signed_corner_areas(corner_areas, coordinates);
 
-  double aspect_ratio = .5 * hm * (a1 + b1 + c1 + d1) / denominator;
+  double aspect_ratio = hm * (a1 + b1 + c1 + d1) / (corner_areas[0] + corner_areas[1] + corner_areas[2] + corner_areas[3]);
 
   if (aspect_ratio > 0)
   {
@@ -744,18 +739,103 @@ double quad_warpage(int /*num_nodes*/, const double coordinates[][3])
 
   jacobian at quad center
  */
-double quad_area(int /*num_nodes*/, const double coordinates[][3])
+double quad_area(int num_nodes, const double coordinates[][3])
 {
-  double corner_areas[4];
-  signed_corner_areas(corner_areas, coordinates);
+  if (4 == num_nodes)
+  {
+    double corner_areas[4];
+    signed_corner_areas(corner_areas, coordinates);
 
-  double area = 0.25 * (corner_areas[0] + corner_areas[1] + corner_areas[2] + corner_areas[3]);
+    double area = 0.25 * (corner_areas[0] + corner_areas[1] + corner_areas[2] + corner_areas[3]);
 
-  if (area > 0)
+    if (area > 0)
+    {
+      return (double)std::min(area, VERDICT_DBL_MAX);
+    }
+    return (double)std::max(area, -VERDICT_DBL_MAX);
+  }
+  else 
   {
-    return (double)std::min(area, VERDICT_DBL_MAX);
+    double area = 0;
+    double tmp_coords[4][3];
+    
+    if (5 == num_nodes)
+    {
+      std::vector< std::vector<int> > tri_conn = { {0,1}, {1,2}, {2,3}, {3,0} };
+
+      //center node 4
+      tmp_coords[2][0] = coordinates[4][0];
+      tmp_coords[2][1] = coordinates[4][1];
+      tmp_coords[2][2] = coordinates[4][2];
+
+      for (std::vector<int> v : tri_conn)
+      {
+        tmp_coords[0][0] = coordinates[v[0]][0];
+        tmp_coords[0][1] = coordinates[v[0]][1];
+        tmp_coords[0][2] = coordinates[v[0]][2];
+        tmp_coords[1][0] = coordinates[v[1]][0];
+        tmp_coords[1][1] = coordinates[v[1]][1];
+        tmp_coords[1][2] = coordinates[v[1]][2];
+        area += tri_area(3, tmp_coords);
+      }     
+    }
+    else if (8 == num_nodes)
+    {
+      std::vector< std::vector<int> > tri_conn = 
+        { {0,4,7}, {4,1,5}, {5,2,6}, {6,3,7} };
+
+      for (std::vector<int> v : tri_conn)
+      {
+        tmp_coords[0][0] = coordinates[v[0]][0];
+        tmp_coords[0][1] = coordinates[v[0]][1];
+        tmp_coords[0][2] = coordinates[v[0]][2];
+        tmp_coords[1][0] = coordinates[v[1]][0];
+        tmp_coords[1][1] = coordinates[v[1]][1];
+        tmp_coords[1][2] = coordinates[v[1]][2];
+        tmp_coords[2][0] = coordinates[v[2]][0];
+        tmp_coords[2][1] = coordinates[v[2]][1];
+        tmp_coords[2][2] = coordinates[v[2]][2];
+        area += tri_area(3, tmp_coords);
+      }
+
+      //interior quad 4567
+      tmp_coords[0][0] = coordinates[4][0];
+      tmp_coords[0][1] = coordinates[4][1];
+      tmp_coords[0][2] = coordinates[4][2];
+      tmp_coords[1][0] = coordinates[5][0];
+      tmp_coords[1][1] = coordinates[5][1];
+      tmp_coords[1][2] = coordinates[5][2];
+      tmp_coords[2][0] = coordinates[6][0];
+      tmp_coords[2][1] = coordinates[6][1];
+      tmp_coords[2][2] = coordinates[6][2];
+      tmp_coords[3][0] = coordinates[7][0];
+      tmp_coords[3][1] = coordinates[7][1];
+      tmp_coords[3][2] = coordinates[7][2];
+      area += quad_area(4, tmp_coords);
+    }
+    else if (9 == num_nodes)
+    {
+      std::vector< std::vector<int> > tri_conn =
+      { {0,4}, {4,1}, {1,5}, {5,2}, {2,6}, {6,3}, {3,7}, {7,0} };
+
+      //quad center node
+      tmp_coords[2][0] = coordinates[8][0];
+      tmp_coords[2][1] = coordinates[8][1];
+      tmp_coords[2][2] = coordinates[8][2];
+
+      for (std::vector<int> v : tri_conn)
+      {
+        tmp_coords[0][0] = coordinates[v[0]][0];
+        tmp_coords[0][1] = coordinates[v[0]][1];
+        tmp_coords[0][2] = coordinates[v[0]][2];
+        tmp_coords[1][0] = coordinates[v[1]][0];
+        tmp_coords[1][1] = coordinates[v[1]][1];
+        tmp_coords[1][2] = coordinates[v[1]][2];
+        area += tri_area(3, tmp_coords);
+      }
+    }
+    return area;
   }
-  return (double)std::max(area, -VERDICT_DBL_MAX);
 }
 
 /*!

ThirdParty/verdict/vtkverdict/V_TetMetric.cpp

@@ -278,20 +278,16 @@ double tet_edge_ratio(int /*num_nodes*/, const double coordinates[][3])
   minimum of the jacobian divided by the lengths of 3 edge vectors
 
  */
-double tet_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
+
+template <typename CoordsContainerType>
+double tet_scaled_jacobian_impl(int /*num_nodes*/, const CoordsContainerType coordinates)
 {
-  const VerdictVector side0(coordinates[1][0] - coordinates[0][0],
-    coordinates[1][1] - coordinates[0][1], coordinates[1][2] - coordinates[0][2]);
-  const VerdictVector side1(coordinates[2][0] - coordinates[1][0],
-    coordinates[2][1] - coordinates[1][1], coordinates[2][2] - coordinates[1][2]);
-  const VerdictVector side2(coordinates[0][0] - coordinates[2][0],
-    coordinates[0][1] - coordinates[2][1], coordinates[0][2] - coordinates[2][2]);
-  const VerdictVector side3(coordinates[3][0] - coordinates[0][0],
-    coordinates[3][1] - coordinates[0][1], coordinates[3][2] - coordinates[0][2]);
-  const VerdictVector side4(coordinates[3][0] - coordinates[1][0],
-    coordinates[3][1] - coordinates[1][1], coordinates[3][2] - coordinates[1][2]);
-  const VerdictVector side5(coordinates[3][0] - coordinates[2][0],
-    coordinates[3][1] - coordinates[2][1], coordinates[3][2] - coordinates[2][2]);
+  const VerdictVector side0{coordinates[0], coordinates[1]};
+  const VerdictVector side1{coordinates[1], coordinates[2]};
+  const VerdictVector side2{coordinates[2], coordinates[0]};
+  const VerdictVector side3{coordinates[0], coordinates[3]};
+  const VerdictVector side4{coordinates[1], coordinates[3]};
+  const VerdictVector side5{coordinates[2], coordinates[3]};
 
   const double jacobi = side3 % (side2 * side0);
 
@@ -303,8 +299,8 @@ double tet_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
   const double side4_length_squared = side4.length_squared();
   const double side5_length_squared = side5.length_squared();
 
-  const double length_squared[4] = { side0_length_squared * side2_length_squared *
-      side3_length_squared,
+  const double length_squared[4] = {
+    side0_length_squared * side2_length_squared * side3_length_squared,
     side0_length_squared * side1_length_squared * side4_length_squared,
     side1_length_squared * side2_length_squared * side5_length_squared,
     side3_length_squared * side4_length_squared * side5_length_squared };
@@ -336,6 +332,15 @@ double tet_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
   return (double)(sqrt2 * jacobi / length_product);
 }
 
+double tet_scaled_jacobian(int num_nodes, const double coordinates[][3])
+{
+    return tet_scaled_jacobian_impl(num_nodes, coordinates);
+}
+
+double tet_scaled_jacobian_from_loc_ptrs(int num_nodes, const double * const * coordinates)
+{
+    return tet_scaled_jacobian_impl(num_nodes, coordinates);
+}
 /*!
   The radius ratio of a tet
 
@@ -401,19 +406,13 @@ double tet_radius_ratio(int /*num_nodes*/, const double coordinates[][3])
     conditioned. Previously, this would only happen when the volume was small
     and positive, but now ill-conditioned inverted tetrahedra are also included.
  */
-double tet_aspect_ratio(int /*num_nodes*/, const double coordinates[][3])
+template <typename CoordsContainerType>
+double tet_aspect_ratio_impl(int /*num_nodes*/, const CoordsContainerType coordinates)
 {
   // Determine side vectors
-  VerdictVector ab, bc, ac, ad, bd, cd;
-
-  ab.set(coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
-    coordinates[1][2] - coordinates[0][2]);
-
-  ac.set(coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
-    coordinates[2][2] - coordinates[0][2]);
-
-  ad.set(coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
-    coordinates[3][2] - coordinates[0][2]);
+  const VerdictVector ab{coordinates[0], coordinates[1]};
+  const VerdictVector ac{coordinates[0], coordinates[2]};
+  const VerdictVector ad{coordinates[0], coordinates[3]};
 
   double detTet = ab % (ac * ad);
 
@@ -422,27 +421,22 @@ double tet_aspect_ratio(int /*num_nodes*/, const double coordinates[][3])
     return (double)VERDICT_DBL_MAX;
   }
 
-  bc.set(coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
-    coordinates[2][2] - coordinates[1][2]);
-
-  bd.set(coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
-    coordinates[3][2] - coordinates[1][2]);
+  VerdictVector bc{coordinates[1], coordinates[2]};
+  VerdictVector bd{coordinates[1], coordinates[3]};
+  VerdictVector cd{coordinates[2], coordinates[3]};
 
-  cd.set(coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
-    coordinates[3][2] - coordinates[2][2]);
-
-  double ab2 = ab.length_squared();
-  double bc2 = bc.length_squared();
-  double ac2 = ac.length_squared();
-  double ad2 = ad.length_squared();
-  double bd2 = bd.length_squared();
-  double cd2 = cd.length_squared();
+  const double ab2 = ab.length_squared();
+  const double bc2 = bc.length_squared();
+  const double ac2 = ac.length_squared();
+  const double ad2 = ad.length_squared();
+  const double bd2 = bd.length_squared();
+  const double cd2 = cd.length_squared();
 
   double A = ab2 > bc2 ? ab2 : bc2;
   double B = ac2 > ad2 ? ac2 : ad2;
   double C = bd2 > cd2 ? bd2 : cd2;
   double D = A > B ? A : B;
-  double hm = D > C ? std::sqrt(D) : std::sqrt(C);
+  const double hm = D > C ? std::sqrt(D) : std::sqrt(C);
 
   bd = ab * bc;
   A = bd.length();
@@ -458,6 +452,15 @@ double tet_aspect_ratio(int /*num_nodes*/, const double coordinates[][3])
   return fix_range(aspect_ratio);
 }
 
+double tet_aspect_ratio(int num_nodes, const double coordinates[][3])
+{
+    return tet_aspect_ratio_impl(num_nodes, coordinates);
+}
+double tet_aspect_ratio_from_loc_ptrs(int num_nodes, const double * const *coordinates)
+{
+    return tet_aspect_ratio_impl(num_nodes, coordinates);
+}
+
 /*!
   the aspect gamma of a tet
 
@@ -954,25 +957,21 @@ double calculate_tet_volume_using_sides(
 
   1/6 * jacobian at a corner node
  */
-double tet_volume(int num_nodes, const double coordinates[][3])
+template <typename CoordsContainerType>
+double tet_volume_impl(int num_nodes, const CoordsContainerType coordinates)
 {
   // Determine side vectors
-  VerdictVector side0, side2, side3;
   if (4 == num_nodes)
   {
-    side2.set(coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
-      coordinates[1][2] - coordinates[0][2]);
-
-    side0.set(coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
-      coordinates[2][2] - coordinates[0][2]);
-
-    side3.set(coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
-      coordinates[3][2] - coordinates[0][2]);
+    const VerdictVector side2{coordinates[0], coordinates[1]};
+    const VerdictVector side0{coordinates[0], coordinates[2]};
+    const VerdictVector side3{coordinates[0], coordinates[3]};
     return calculate_tet_volume_using_sides(side0, side2, side3);
   }
   else
   {
     VerdictVector tet_pts[15];
+    VerdictVector side0, side2, side3;
 
     // create a vector for each point
     for (int k = 0; k < num_nodes; k++)
@@ -981,7 +980,7 @@ double tet_volume(int num_nodes, const double coordinates[][3])
     }
 
     // determine center point of the higher-order nodes
-    VerdictVector centroid(0, 0, 0);
+    VerdictVector centroid(0.0, 0.0, 0.0);
     for (int k = 4; k < num_nodes; k++)
     {
       centroid += VerdictVector(coordinates[k][0], coordinates[k][1], coordinates[k][2]);
@@ -1133,6 +1132,16 @@ double tet_volume(int num_nodes, const double coordinates[][3])
   return 0;
 }
 
+double tet_volume(int num_nodes, const double coordinates[][3])
+{
+    return tet_volume_impl(num_nodes, coordinates);
+}
+
+double tet_volume_from_loc_ptrs(int num_nodes, const double * const *coordinates)
+{
+    return tet_volume_impl(num_nodes, coordinates);
+}
+
 /*!
   the condition of a tet
 
@@ -1143,30 +1152,20 @@ double tet_volume(int num_nodes, const double coordinates[][3])
     conditioned. Previously, this would only happen when the volume was small
     and positive, but now ill-conditioned inverted tetrahedra are also included.
  */
-double tet_condition(int /*num_nodes*/, const double coordinates[][3])
+template <typename CoordsContainerType>
+double tet_condition_impl(int /*num_nodes*/, const CoordsContainerType coordinates)
 {
-  double condition, term1, term2, det;
-
-  VerdictVector side0, side2, side3;
-
-  side0.set(coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
-    coordinates[1][2] - coordinates[0][2]);
+  const VerdictVector side0{coordinates[0], coordinates[1]};
+  const VerdictVector side2{coordinates[2], coordinates[0]};
+  const VerdictVector side3{coordinates[0], coordinates[3]};
 
-  side2.set(coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
-    coordinates[0][2] - coordinates[2][2]);
-
-  side3.set(coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
-    coordinates[3][2] - coordinates[0][2]);
-
-  VerdictVector c_1, c_2, c_3;
+  const VerdictVector c_1 = side0;
+  const VerdictVector c_2 = (-2 * side2 - side0) / sqrt3;
+  const VerdictVector c_3 = (3 * side3 + side2 - side0) / sqrt6;
 
-  c_1 = side0;
-  c_2 = (-2 * side2 - side0) / sqrt3;
-  c_3 = (3 * side3 + side2 - side0) / sqrt6;
-
-  term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
-  term2 = (c_1 * c_2) % (c_1 * c_2) + (c_2 * c_3) % (c_2 * c_3) + (c_1 * c_3) % (c_1 * c_3);
-  det = c_1 % (c_2 * c_3);
+  const double term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
+  const double term2 = (c_1 * c_2) % (c_1 * c_2) + (c_2 * c_3) % (c_2 * c_3) + (c_1 * c_3) % (c_1 * c_3);
+  const double det = c_1 % (c_2 * c_3);
 
   if (std::abs(det) <= VERDICT_DBL_MIN)
   {
@@ -1174,12 +1173,19 @@ double tet_condition(int /*num_nodes*/, const double coordinates[][3])
   }
   else
   {
-    condition = std::sqrt(term1 * term2) / (3.0 * det);
+    return std::sqrt(term1 * term2) / (3.0 * det);
   }
+}
 
-  return (double)condition;
+double tet_condition(int num_nodes, const double coordinates[][3])
+{
+    return tet_condition_impl(num_nodes, coordinates);
 }
 
+double tet_condition_from_loc_ptrs(int num_nodes, const double * const *coordinates)
+{
+    return tet_condition_impl(num_nodes, coordinates);
+}
 /*!
   the jacobian of a tet
 
@@ -1517,9 +1523,10 @@ double tet_timestep(int num_nodes, const double coordinates[][3], double density
   return char_length / denominator;
 }
 
-VerdictVector tet10_auxillary_node_coordinate(const double coordinates[][3])
+template <typename CoordsContainerType>
+VerdictVector tet10_auxillary_node_coordinate(const CoordsContainerType coordinates)
 {
-  VerdictVector aux_node(0, 0, 0);
+  VerdictVector aux_node(0.0, 0.0, 0.0);
   for (int i = 4; i < 10; i++)
   {
     VerdictVector tmp_vec(coordinates[i][0], coordinates[i][1], coordinates[i][2]);
@@ -1530,7 +1537,8 @@ VerdictVector tet10_auxillary_node_coordinate(const double coordinates[][3])
   return aux_node;
 }
 
-double tet10_min_inradius(const double coordinates[][3], int begin_index, int end_index)
+template <typename CoordsContainerType>
+double tet10_min_inradius(const CoordsContainerType coordinates, int begin_index, int end_index)
 {
   double min_tetinradius = VERDICT_DBL_MAX;
 
@@ -1582,7 +1590,9 @@ double tet10_characteristic_length(const double coordinates[][3])
 
   return min_tetinradius;
 }
-double calculate_tet4_outer_radius(const double coordinates[][3])
+
+template <typename CoordsContainerType>
+double calculate_tet4_outer_radius(const CoordsContainerType coordinates)
 {
   verdict::VerdictVector nE[4];
   for (int i{ 0 }; i < 4; i++)
@@ -1604,7 +1614,8 @@ double calculate_tet4_outer_radius(const double coordinates[][3])
   return outer_radius;
 }
 
-double tet10_normalized_inradius(const double coordinates[][3])
+template <typename CoordsContainerType>
+double tet10_normalized_inradius(const CoordsContainerType coordinates)
 {
   double min_inradius_for_subtet_with_parent_node = tet10_min_inradius(coordinates, 0, 3);
   double min_inradius_for_subtet_with_no_parent_node = tet10_min_inradius(coordinates, 4, 11);
@@ -1621,7 +1632,8 @@ double tet10_normalized_inradius(const double coordinates[][3])
   return fix_range(norm_inrad);
 }
 
-double tet4_normalized_inradius(const double coordinates[][3])
+template <typename CoordsContainerType>
+double tet4_normalized_inradius(const CoordsContainerType coordinates)
 {
   double tet10_coords[10][3];
   for (int i = 0; i < 4; i++)
@@ -1643,7 +1655,8 @@ double tet4_normalized_inradius(const double coordinates[][3])
   return tet10_normalized_inradius(tet10_coords);
 }
 
-double tet_normalized_inradius(int num_nodes, const double coordinates[][3])
+template <typename CoordsContainerType>
+double tet_normalized_inradius_impl(int num_nodes, const CoordsContainerType coordinates)
 {
   if (num_nodes == 4)
   {
@@ -1656,31 +1669,32 @@ double tet_normalized_inradius(int num_nodes, const double coordinates[][3])
   return 0.0;
 }
 
-double tet_mean_ratio(int /*num_nodes*/, const double coordinates[][3])
+double tet_normalized_inradius(int num_nodes, const double coordinates[][3])
 {
-  const VerdictVector side0(coordinates[1][0] - coordinates[0][0],
-    coordinates[1][1] - coordinates[0][1], coordinates[1][2] - coordinates[0][2]);
+    return tet_normalized_inradius_impl(num_nodes, coordinates);
+}
 
-  const VerdictVector side2(coordinates[0][0] - coordinates[2][0],
-    coordinates[0][1] - coordinates[2][1], coordinates[0][2] - coordinates[2][2]);
+double tet_normalized_inradius_from_loc_ptrs(int num_nodes, const double * const *coordinates)
+{
+    return tet_normalized_inradius_impl(num_nodes, coordinates);
+}
 
-  const VerdictVector side3(coordinates[3][0] - coordinates[0][0],
-    coordinates[3][1] - coordinates[0][1], coordinates[3][2] - coordinates[0][2]);
+template <typename CoordsContainerType>
+double tet4_mean_ratio(const CoordsContainerType coordinates)
+{
+  const VerdictVector side0{coordinates[0], coordinates[1]};
+  const VerdictVector side2{coordinates[2], coordinates[0]};
+  const VerdictVector side3{coordinates[0], coordinates[3]};
 
   const double tetVolume = calculate_tet_volume_using_sides(side0, side2, side3);
-  if (std::abs(tetVolume) < VERDICT_DBL_MIN)
+  if (std::abs( tetVolume ) < VERDICT_DBL_MIN)
   {
     return 0.0;
   }
 
-  const VerdictVector side1(coordinates[2][0] - coordinates[1][0],
-    coordinates[2][1] - coordinates[1][1], coordinates[2][2] - coordinates[1][2]);
-
-  const VerdictVector side4(coordinates[3][0] - coordinates[1][0],
-    coordinates[3][1] - coordinates[1][1], coordinates[3][2] - coordinates[1][2]);
-
-  const VerdictVector side5(coordinates[3][0] - coordinates[2][0],
-    coordinates[3][1] - coordinates[2][1], coordinates[3][2] - coordinates[2][2]);
+  const VerdictVector side1{coordinates[1], coordinates[2]};
+  const VerdictVector side4{coordinates[1], coordinates[3]};
+  const VerdictVector side5{coordinates[2], coordinates[3]};
 
   const double side0_length_squared = side0.length_squared();
   const double side1_length_squared = side1.length_squared();
@@ -1689,9 +1703,86 @@ double tet_mean_ratio(int /*num_nodes*/, const double coordinates[][3])
   const double side4_length_squared = side4.length_squared();
   const double side5_length_squared = side5.length_squared();
 
-  const int sign = tetVolume < 0. ? -1 : 1;
-  return sign * 12. * std::pow(3. * std::abs(tetVolume), 2. / 3.) /
-    (side0_length_squared + side1_length_squared + side2_length_squared + side3_length_squared +
-      side4_length_squared + side5_length_squared);
+  //const int sign = tetVolume < 0. ? -1 : 1;
+  //return sign * 12. * std::pow(3.*fabs(tetVolume), 2./3.) / (side0_length_squared + side1_length_squared + side2_length_squared + side3_length_squared + side4_length_squared + side5_length_squared);
+  double sum = (side0_length_squared + side1_length_squared + side2_length_squared +
+    side3_length_squared + side4_length_squared + side5_length_squared) / 6;
+  return 6 * std::pow(2, 0.5) * tetVolume / std::pow(sum, 3. / 2.);    
+}
+
+template <typename CoordsContainerType>
+double tet10_mean_ratio(const CoordsContainerType coordinates)
+{
+  double min_tet_mean_ratio = VERDICT_DBL_MAX;
+
+  VerdictVector auxillary_node = tet10_auxillary_node_coordinate(coordinates);
+
+  double aux_node_scale = 3.0 * std::pow(3., 0.5) * 0.25;
+
+  for (int i = 0; i <= 11; i++)
+  {
+    int subtet_conn[4];
+    subtet_conn[0] = tet10_subtet_conn[i][0];
+    subtet_conn[1] = tet10_subtet_conn[i][1];
+    subtet_conn[2] = tet10_subtet_conn[i][2];
+    subtet_conn[3] = tet10_subtet_conn[i][3];
+
+    //get the coordinates of the nodes
+    double subtet_coords[4][3];
+    for (int k = 0; k < 4; k++)
+    {
+      int node_index = subtet_conn[k];
+
+      if (10 == node_index)
+      {
+        subtet_coords[k][0] = auxillary_node.x();
+        subtet_coords[k][1] = auxillary_node.y();
+        subtet_coords[k][2] = auxillary_node.z();
+      }
+      else
+      {
+        subtet_coords[k][0] = coordinates[node_index][0];
+        subtet_coords[k][1] = coordinates[node_index][1];
+        subtet_coords[k][2] = coordinates[node_index][2];
+      }
+    }
+
+    double tmp_mean_ratio = tet4_mean_ratio(subtet_coords);
+
+    if (i > 3)
+    {
+      tmp_mean_ratio *= aux_node_scale;
+    }
+
+    if (tmp_mean_ratio < min_tet_mean_ratio)
+    {
+      min_tet_mean_ratio = tmp_mean_ratio;
+    }
+  }
+  return min_tet_mean_ratio;
+}
+
+template <typename CoordsContainerType>
+double tet_mean_ratio_impl(int num_nodes, const CoordsContainerType coordinates)
+{
+  if (num_nodes == 4)
+  {
+    return tet4_mean_ratio(coordinates);
+  }
+  else if (num_nodes >= 10)
+  {
+    return tet10_mean_ratio(coordinates);
+  }
+  return 0.0;
+}
+
+double tet_mean_ratio(int num_nodes, const double coordinates[][3])
+{
+   return tet_mean_ratio_impl(num_nodes, coordinates);
+}
+
+double tet_mean_ratio_from_loc_ptrs(int num_nodes, const double * const * coordinates)
+{
+   return tet_mean_ratio_impl(num_nodes, coordinates);
 }
 } // namespace verdict

ThirdParty/verdict/vtkverdict/V_TriMetric.cpp

@@ -24,6 +24,7 @@
 #include "verdict_defines.hpp"
 
 #include <algorithm>
+#include <array>
 
 namespace VERDICT_NAMESPACE
 {
@@ -149,27 +150,23 @@ double tri_edge_ratio(int /*num_nodes*/, const double coordinates[][3])
      what is now called "v_tri_aspect_frobenius"
 
  */
-double tri_aspect_ratio(int /*num_nodes*/, const double coordinates[][3])
+template <typename CoordsContainerType>
+double tri_aspect_ratio_impl(int /*num_nodes*/, const CoordsContainerType coordinates, const int dimension)
 {
   // three vectors for each side
-  VerdictVector a(coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
-    coordinates[1][2] - coordinates[0][2]);
-
-  VerdictVector b(coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
-    coordinates[2][2] - coordinates[1][2]);
+  const VerdictVector a{coordinates[0], coordinates[1], dimension};
+  const VerdictVector b{coordinates[1], coordinates[2], dimension};
+  const VerdictVector c{coordinates[2], coordinates[0], dimension};
 
-  VerdictVector c(coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
-    coordinates[0][2] - coordinates[2][2]);
-
-  double a1 = a.length();
-  double b1 = b.length();
-  double c1 = c.length();
+  const double a1 = a.length();
+  const double b1 = b.length();
+  const double c1 = c.length();
 
   double hm = a1 > b1 ? a1 : b1;
   hm = hm > c1 ? hm : c1;
 
-  VerdictVector ab = a * b;
-  double denominator = ab.length();
+  const VerdictVector ab = a * b;
+  const double denominator = ab.length();
 
   if (denominator < VERDICT_DBL_MIN)
   {
@@ -177,9 +174,7 @@ double tri_aspect_ratio(int /*num_nodes*/, const double coordinates[][3])
   }
   else
   {
-    double aspect_ratio;
-    aspect_ratio = aspect_ratio_normal_coeff * hm * (a1 + b1 + c1) / denominator;
-
+    const double aspect_ratio = aspect_ratio_normal_coeff * hm * (a1 + b1 + c1) / denominator;
     if (aspect_ratio > 0)
     {
       return (double)std::min(aspect_ratio, VERDICT_DBL_MAX);
@@ -188,6 +183,14 @@ double tri_aspect_ratio(int /*num_nodes*/, const double coordinates[][3])
   }
 }
 
+double tri_aspect_ratio(int num_nodes, const double coordinates[][3])
+{
+    return tri_aspect_ratio_impl(num_nodes, coordinates, 3);
+}
+double tri_aspect_ratio_from_loc_ptrs(int num_nodes, const double * const *coordinates, const int dimension)
+{
+    return tri_aspect_ratio_impl(num_nodes, coordinates, dimension);
+}
 /*!
    the radius ratio of a triangle
 
@@ -279,26 +282,124 @@ double tri_aspect_frobenius(int /*num_nodes*/, const double coordinates[][3])
 
   0.5 * jacobian at a node
  */
-double tri_area(int /*num_nodes*/, const double coordinates[][3])
+template <typename CoordsContainerType>
+double tri_area_impl(int num_nodes, const CoordsContainerType coordinates, const int dimension)
 {
-  // two vectors for two sides
-  VerdictVector side1(coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
-    coordinates[1][2] - coordinates[0][2]);
-
-  VerdictVector side3(coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
-    coordinates[2][2] - coordinates[0][2]);
+  if (3 == num_nodes)
+  {
+    // two vectors for two sides
+    const VerdictVector side1{ coordinates[0], coordinates[1], dimension };
+    const VerdictVector side3{ coordinates[0], coordinates[2], dimension };
 
-  // the cross product of the two vectors representing two sides of the
-  // triangle
-  VerdictVector tmp = side1 * side3;
+    // the cross product of the two vectors representing two sides of the
+    // triangle
+    const VerdictVector tmp = side1 * side3;
 
-  // return the magnitude of the vector divided by two
-  double area = 0.5 * tmp.length();
-  if (area > 0)
+    // return the magnitude of the vector divided by two
+    const double area = 0.5 * tmp.length();
+    if (area > 0)
+    {
+      return (double)std::min(area, VERDICT_DBL_MAX);
+    }
+    return (double)std::max(area, -VERDICT_DBL_MAX);
+  }
+  else
   {
-    return (double)std::min(area, VERDICT_DBL_MAX);
+    const double *tmp_coords[3];
+    double tri_area = 0.0;
+
+    if (6 == num_nodes)
+    {
+      // 035
+      tmp_coords[0] = coordinates[0];
+      tmp_coords[1] = coordinates[3];
+      tmp_coords[2] = coordinates[5];      
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+
+      // 314
+      tmp_coords[0] = coordinates[3];
+      tmp_coords[1] = coordinates[1];
+      tmp_coords[2] = coordinates[4];      
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+
+      // 425
+      tmp_coords[0] = coordinates[4];
+      tmp_coords[1] = coordinates[2];
+      tmp_coords[2] = coordinates[5];
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+
+      // 345
+      tmp_coords[0] = coordinates[3];
+      tmp_coords[1] = coordinates[4];     
+      tmp_coords[2] = coordinates[5];
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+    }    
+    else if (7 == num_nodes)
+    {
+      //center node 6
+      tmp_coords[2] = coordinates[6];
+      
+      // 036
+      tmp_coords[0] = coordinates[0];
+      tmp_coords[1] = coordinates[3];      
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+
+      // 316
+      tmp_coords[0] = coordinates[3];
+      tmp_coords[1] = coordinates[1];
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+
+      // 146
+      tmp_coords[0] = coordinates[1];
+      tmp_coords[1] = coordinates[4];
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+
+      // 426
+      tmp_coords[0] = coordinates[4];
+      tmp_coords[1] = coordinates[2];
+      tri_area += tri_area_impl(3, tmp_coords, dimension);      
+
+      // 256
+      tmp_coords[0] = coordinates[2];
+      tmp_coords[1] = coordinates[5];
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+
+      // 506
+      tmp_coords[0] = coordinates[5];
+      tmp_coords[1] = coordinates[0];
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+    }    
+    else if( 4 == num_nodes )
+    {
+      //center node 3
+      tmp_coords[2] = coordinates[3];      
+
+      // 013
+      tmp_coords[0] = coordinates[0];
+      tmp_coords[1] = coordinates[1];      
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+
+      // 123
+      tmp_coords[0] = coordinates[1];
+      tmp_coords[1] = coordinates[2];
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+
+      // 203
+      tmp_coords[0] = coordinates[2];
+      tmp_coords[1] = coordinates[0];
+      tri_area += tri_area_impl(3, tmp_coords, dimension);
+    }    
+    return tri_area;
   }
-  return (double)std::max(area, -VERDICT_DBL_MAX);
+}
+double tri_area(int num_nodes, const double coordinates[][3])
+{
+    return tri_area_impl(num_nodes, coordinates, 3);
+}
+
+double tri_area_from_loc_ptrs(int num_nodes, const double * const *coordinates, const int dimension)
+{
+    return tri_area_impl(num_nodes, coordinates, dimension);
 }
 
 /*!
@@ -463,54 +564,56 @@ double tri_equiangle_skew(int num_nodes, const double coordinates[][3])
 
   Condition number of the jacobian matrix at any corner
  */
-double tri_condition(int /*num_nodes*/, const double coordinates[][3])
+template <typename CoordsContainerType>
+double tri_condition_impl(int /*num_nodes*/, const CoordsContainerType coordinates, const int dimension)
 {
-  VerdictVector v1(coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
-    coordinates[1][2] - coordinates[0][2]);
-
-  VerdictVector v2(coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
-    coordinates[2][2] - coordinates[0][2]);
-
-  VerdictVector tri_normal = v1 * v2;
-  double areax2 = tri_normal.length();
+  const VerdictVector v1{coordinates[0], coordinates[1], dimension};
+  const VerdictVector v2{coordinates[0], coordinates[2], dimension};
+  const VerdictVector tri_normal = v1 * v2;
+  const double areax2 = tri_normal.length();
 
   if (areax2 == 0.0)
   {
     return (double)VERDICT_DBL_MAX;
   }
 
-  double condition = (double)(((v1 % v1) + (v2 % v2) - (v1 % v2)) / (areax2 * sqrt3));
+  const double condition = (double)(((v1 % v1) + (v2 % v2) - (v1 % v2)) / (areax2 * sqrt3));
 
   return (double)std::min(condition, VERDICT_DBL_MAX);
 }
 
+double tri_condition(int num_nodes, const double coordinates[][3])
+{
+    return tri_condition_impl(num_nodes, coordinates, 3);
+}
+
+double tri_condition_from_loc_ptrs(int num_nodes, const double * const * coordinates, const int dimension)
+{
+    return tri_condition_impl(num_nodes, coordinates, dimension);
+}
+
 /*!
   The scaled jacobian of a tri
 
   minimum of the jacobian divided by the lengths of 2 edge vectors
  */
-double tri_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
+template <typename CoordsContainerType>
+double tri_scaled_jacobian_impl(int /*num_nodes*/, const CoordsContainerType coordinates, const int dimension)
 {
-  VerdictVector first, second;
-  double jacobian;
-
-  VerdictVector edge[3];
-  edge[0].set(coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
-    coordinates[1][2] - coordinates[0][2]);
-
-  edge[1].set(coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
-    coordinates[2][2] - coordinates[0][2]);
+  const VerdictVector edge[3] =
+  {
+      {coordinates[0], coordinates[1], dimension},
+      {coordinates[0], coordinates[2], dimension},
+      {coordinates[1], coordinates[2], dimension},
+  };
 
-  edge[2].set(coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
-    coordinates[2][2] - coordinates[1][2]);
-  first = edge[1] - edge[0];
-  second = edge[2] - edge[0];
+  const VerdictVector first = edge[1] - edge[0];
+  const VerdictVector second = edge[2] - edge[0];
 
-  VerdictVector cross = first * second;
-  jacobian = cross.length();
+  const VerdictVector cross = first * second;
+  double jacobian = cross.length();
 
-  double max_edge_length_product;
-  max_edge_length_product = std::max(edge[0].length() * edge[1].length(),
+  const double max_edge_length_product = std::max(edge[0].length() * edge[1].length(),
     std::max(edge[1].length() * edge[2].length(), edge[0].length() * edge[2].length()));
 
   if (max_edge_length_product < VERDICT_DBL_MIN)
@@ -528,6 +631,16 @@ double tri_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
   return (double)std::max(jacobian, -VERDICT_DBL_MAX);
 }
 
+double tri_scaled_jacobian(int num_nodes, const double coordinates[][3])
+{
+    return tri_scaled_jacobian_impl(num_nodes, coordinates, 3);
+}
+
+double tri_scaled_jacobian_from_loc_ptrs(int num_nodes, const double * const * coordinates, const int dimension)
+{
+    return tri_scaled_jacobian_impl(num_nodes, coordinates, dimension);
+}
+
 /*!
   The shape of a tri
 
@@ -649,10 +762,11 @@ double tri_distortion(int num_nodes, const double coordinates[][3])
     return (double)distortion;
   }
 
-  else if (num_nodes == 6)
+  else if (num_nodes >= 6)
   {
     total_number_of_gauss_points = 6;
     number_of_gauss_points = 6;
+    num_nodes = 6; //seven nodes not handled 
   }
 
   distortion = VERDICT_DBL_MAX;
@@ -829,7 +943,8 @@ double tri_distortion(int num_nodes, const double coordinates[][3])
   return (double)std::max(distortion, -VERDICT_DBL_MAX);
 }
 
-double tri_inradius(const double coordinates[][3])
+template <typename CoordsContainerType>
+double tri_inradius(const CoordsContainerType coordinates)
 {
   double sp = 0.0;
   VerdictVector sides[3];
@@ -849,7 +964,8 @@ double tri_inradius(const double coordinates[][3])
   return ir;
 }
 
-double tri6_min_inradius(const double coordinates[][3])
+template <typename CoordsContainerType>
+double tri6_min_inradius(const CoordsContainerType coordinates, const int dimension)
 {
   static int SUBTRI_NODES[4][3] = { { 0, 3, 5 }, { 3, 1, 4 }, { 5, 4, 2 }, { 3, 4, 5 } };
   double min_inrad = VERDICT_DBL_MAX;
@@ -861,7 +977,7 @@ double tri6_min_inradius(const double coordinates[][3])
       int idx = SUBTRI_NODES[i][j];
       subtri_coords[j][0] = coordinates[idx][0];
       subtri_coords[j][1] = coordinates[idx][1];
-      subtri_coords[j][2] = coordinates[idx][2];
+      subtri_coords[j][2] = dimension == 2 ? 0.0 : coordinates[idx][2];
     }
     double subtri_inrad = tri_inradius(subtri_coords);
     if (subtri_inrad < min_inrad)
@@ -872,7 +988,8 @@ double tri6_min_inradius(const double coordinates[][3])
   return min_inrad;
 }
 
-double calculate_tri3_outer_radius(const double coordinates[][3])
+template <typename CoordsContainerType>
+double calculate_tri3_outer_radius(const CoordsContainerType coordinates, const int dimension)
 {
   double sp = 0.0;
   VerdictVector sides[3];
@@ -880,8 +997,8 @@ double calculate_tri3_outer_radius(const double coordinates[][3])
   for (int i = 0; i < 3; i++)
   {
     int j = (i + 1) % 3;
-    sides[i].set(coordinates[j][0] - coordinates[i][0], coordinates[j][1] - coordinates[i][1],
-      coordinates[j][2] - coordinates[i][2]);
+    const VerdictVector thisSide{coordinates[i], coordinates[j], dimension};
+    sides[i] = thisSide;
     slen[i] = sides[i].length();
     sp += slen[i];
   }
@@ -895,50 +1012,80 @@ double calculate_tri3_outer_radius(const double coordinates[][3])
   return cr;
 }
 
-double tri6_normalized_inradius(const double coordinates[][3])
+static double fix_range(double v)
+{
+  if (std::isnan(v))
+  {
+    return VERDICT_DBL_MAX;
+  }
+  if (v >= VERDICT_DBL_MAX)
+  {
+    return VERDICT_DBL_MAX;
+  }
+  if (v <= -VERDICT_DBL_MAX)
+  {
+    return -VERDICT_DBL_MAX;
+  }
+  return v;
+}
+
+template <typename CoordsContainerType>
+double tri6_normalized_inradius(const CoordsContainerType coordinates, const int dimension)
 {
-  double min_inradius_for_subtri = tri6_min_inradius(coordinates);
-  double outer_radius = calculate_tri3_outer_radius(coordinates);
+  double min_inradius_for_subtri = tri6_min_inradius(coordinates, dimension);
+  double outer_radius = calculate_tri3_outer_radius(coordinates, dimension);
   double normalized_inradius = 4.0 * min_inradius_for_subtri / outer_radius;
 
-  return normalized_inradius;
+  return fix_range(normalized_inradius);
 }
 
-double tri3_normalized_inradius(const double coordinates[][3])
+template <typename CoordsContainerType>
+double tri3_normalized_inradius(const CoordsContainerType coordinates, const int dimension)
 {
-  double tri6_coords[6][3];
+  std::array<const double*, 6> tri6_coords;
   for (int i = 0; i < 3; i++)
-  {
-    tri6_coords[i][0] = coordinates[i][0];
-    tri6_coords[i][1] = coordinates[i][1];
-    tri6_coords[i][2] = coordinates[i][2];
-  }
+    tri6_coords[i] = coordinates[i];
 
   static int eidx[3][2] = { { 0, 1 }, { 1, 2 }, { 2, 0 } };
+  double tri6_midnode_coords[3][3];
   for (int i = 3; i < 6; i++)
   {
     int i0 = eidx[i - 3][0];
     int i1 = eidx[i - 3][1];
-    tri6_coords[i][0] = (coordinates[i0][0] + coordinates[i1][0]) * 0.5;
-    tri6_coords[i][1] = (coordinates[i0][1] + coordinates[i1][1]) * 0.5;
-    tri6_coords[i][2] = (coordinates[i0][2] + coordinates[i1][2]) * 0.5;
+
+    tri6_midnode_coords[i-3][0] = (coordinates[i0][0] + coordinates[i1][0]) * 0.5;
+    tri6_midnode_coords[i-3][1] = (coordinates[i0][1] + coordinates[i1][1]) * 0.5;
+    tri6_midnode_coords[i-3][2] = (coordinates[i0][2] + coordinates[i1][2]) * 0.5;
+    tri6_coords[i] = tri6_midnode_coords[i-3];
   }
-  return tri6_normalized_inradius(tri6_coords);
+  return tri6_normalized_inradius(tri6_coords.data(), dimension);
 }
 
 //! Calculates the minimum normalized inner radius of a tri6
 /** Ratio of the minimum subtri inner radius to tri outer radius*/
 /* Currently, supports tri 6 and 3.*/
-double tri_normalized_inradius(int num_nodes, const double coordinates[][3])
+template <typename CoordsContainerType>
+double tri_normalized_inradius_impl(int num_nodes, const CoordsContainerType coordinates, const int dimension)
 {
   if (num_nodes == 3)
   {
-    return tri3_normalized_inradius(coordinates);
+    return tri3_normalized_inradius(coordinates, dimension);
   }
   else if (num_nodes == 6)
   {
-    return tri6_normalized_inradius(coordinates);
+    return tri6_normalized_inradius(coordinates, dimension);
   }
   return 0.0;
 }
+
+double tri_normalized_inradius(int num_nodes, const double coordinates[][3])
+{
+    return tri_normalized_inradius_impl(num_nodes, coordinates, 3);
+}
+
+double tri_normalized_inradius_from_loc_ptrs(int num_nodes, const double * const *coordinates, const int dimension)
+{
+    return tri_normalized_inradius_impl(num_nodes, coordinates, dimension);
+}
+
 } // namespace verdict

ThirdParty/verdict/vtkverdict/VerdictVector.hpp

@@ -45,9 +45,14 @@ public:
   //- Constructor: create vector from tuple
 
   VerdictVector(const VerdictVector& tail, const VerdictVector& head);
+  VerdictVector(const double *tail, const double *head, int dimension);
+  VerdictVector(const double *tail, const double *head);
   //- Constructor for a VerdictVector starting at tail and pointing
   //- to head.
 
+  template <typename ARG1, typename ARG2, typename ARG3> VerdictVector(ARG1, ARG2, ARG3) = delete;
+  //- define this template to avoid ambiguity between the (double, double, double) and (double *, double *, int) constructors
+
   VerdictVector(const VerdictVector& copy_from); //- Copy Constructor
 
   //- Heading: Set and Inquire Functions
@@ -267,6 +272,20 @@ inline VerdictVector::VerdictVector()
 {
 }
 
+inline VerdictVector::VerdictVector(const double *tail, const double *head, int dimension)
+  : xVal{head[0] - tail[0]}
+  , yVal{head[1] - tail[1]}
+  , zVal{dimension == 2 ? 0.0 : head[2] - tail[2]}
+{
+}
+
+inline VerdictVector::VerdictVector(const double *tail, const double *head)
+  : xVal{head[0] - tail[0]}
+  , yVal{head[1] - tail[1]}
+  , zVal{head[2] - tail[2]}
+{
+}
+
 inline VerdictVector::VerdictVector(const VerdictVector& tail, const VerdictVector& head)
   : xVal(head.xVal - tail.xVal)
   , yVal(head.yVal - tail.yVal)

ThirdParty/verdict/vtkverdict/verdict.h

@@ -286,6 +286,7 @@ VERDICT_EXPORT double tet_radius_ratio(int num_nodes, const double coordinates[]
    length and the inradius of the tetrahedron
    Reference ---  P. Frey and P.-L. George, Meshing, Hermes (2000). */
 VERDICT_EXPORT double tet_aspect_ratio(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tet_aspect_ratio_from_loc_ptrs(int num_nodes, const double * const *coordinates);
 
 //! Calculates tet aspect gamma metric.
 /**  Srms**3 / (8.479670*V) where Srms = sqrt(Sum(Si**2)/6), Si = edge length.
@@ -313,6 +314,7 @@ VERDICT_EXPORT double tet_collapse_ratio(int num_nodes, const double coordinates
    Reference ---  V. N. Parthasarathy et al, A comparison of tetrahedron
    quality measures, Finite Elem. Anal. Des., Vol 15(1993), 255-261. */
 VERDICT_EXPORT double tet_volume(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tet_volume_from_loc_ptrs(int num_nodes, const double * const *coordinates);
 
 //! Calculates tet condition metric.
 /** Condition number of the Jacobian matrix at any corner.
@@ -320,6 +322,7 @@ VERDICT_EXPORT double tet_volume(int num_nodes, const double coordinates[][3]);
    Optimization of the Jacobian Matrix Norm and Associated Quantities,
    Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185. */
 VERDICT_EXPORT double tet_condition(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tet_condition_from_loc_ptrs(int num_nodes, const double * const *coordinates);
 
 //! Calculates tet jacobian.
 /** Minimum pointwise volume at any corner.
@@ -334,6 +337,7 @@ VERDICT_EXPORT double tet_jacobian(int num_nodes, const double coordinates[][3])
    Optimization of the Jacobian Matrix Norm and Associated Quantities,
    Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185. */
 VERDICT_EXPORT double tet_scaled_jacobian(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tet_scaled_jacobian_from_loc_ptrs(int num_nodes, const double * const * coordinates);
 
 //! Calculates tet mean ratio.
 /** Ratio of tet volume to volume of an equilateral tet with the same RMS edge length
@@ -341,11 +345,13 @@ VERDICT_EXPORT double tet_scaled_jacobian(int num_nodes, const double coordinate
    meshes, IJNME 2010 82:843-867 Reference 2 --- Danial Ibanez - PhD Thesis, Conformal Mesh
    Adaptation on Heterogeneous Supercomputers */
 VERDICT_EXPORT double tet_mean_ratio(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tet_mean_ratio_from_loc_ptrs(int num_nodes, const double * const *coordinates);
 
 //! Calculates the minimum normalized inner radius of a tet
 /** Ratio of the minimum subtet inner radius to tet outer radius*/
 /* Currently supports tetra 10 and 4.*/
 VERDICT_EXPORT double tet_normalized_inradius(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tet_normalized_inradius_from_loc_ptrs(int num_nodes, const double * const *coordinates);
 
 //! Calculates tet shape metric.
 /** 3/Mean Ratio of weighted Jacobian matrix.
@@ -618,6 +624,7 @@ VERDICT_EXPORT double tri_edge_ratio(int num_nodes, const double coordinates[][3
     Reference --- P. P. Pebay & T. J. Baker, Analysis of Triangle Quality
     Measures, AMS Math. Comp., 2003, 72(244):1817-1839 */
 VERDICT_EXPORT double tri_aspect_ratio(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tri_aspect_ratio_from_loc_ptrs(int num_nodes, const double * const * coordinates, const int dimension = 3);
 
 //! Calculates triangle metric.
 /** radius ratio
@@ -632,6 +639,7 @@ VERDICT_EXPORT double tri_aspect_frobenius(int num_nodes, const double coordinat
 //! Calculates triangle metric.
 /** Maximum included angle in triangle */
 VERDICT_EXPORT double tri_area(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tri_area_from_loc_ptrs(int num_nodes, const double * const *coordinates, const int dimension = 3);
 
 //! Calculates triangle metric.
 /** Minimum included angle in triangle */
@@ -647,6 +655,7 @@ VERDICT_EXPORT double tri_maximum_angle(int num_nodes, const double coordinates[
    Optimization of the Jacobian Matrix Norm and Associated Quantities,
    Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185. */
 VERDICT_EXPORT double tri_condition(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tri_condition_from_loc_ptrs(int num_nodes, const double * const *coordinates, const int dimension = 3);
 
 //! Calculates triangle metric.
 /** Minimum Jacobian divided by the lengths of 2 edge vectors.
@@ -654,6 +663,7 @@ VERDICT_EXPORT double tri_condition(int num_nodes, const double coordinates[][3]
    Optimization of the Jacobian Matrix Norm and Associated Quantities,
    Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185. */
 VERDICT_EXPORT double tri_scaled_jacobian(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tri_scaled_jacobian_from_loc_ptrs(int num_nodes, const double * const *coordinates, const int dimension=3);
 
 //! Calculates triangle metric.
 /** Min( J, 1/J ), where J is determinant of weighted Jacobian matrix.
@@ -687,6 +697,7 @@ VERDICT_EXPORT double tri_equiangle_skew(int num_nodes, const double coordinates
 /** Ratio of the minimum subtet inner radius to tet outer radius*/
 /* Currently supports tri 6 and 3.*/
 VERDICT_EXPORT double tri_normalized_inradius(int num_nodes, const double coordinates[][3]);
+VERDICT_EXPORT double tri_normalized_inradius_from_loc_ptrs(int num_nodes, const double * const *coordinates, const int dimension=3);
 } // namespace verdict
 
 #endif /* __verdict_h */