Add circumcenter calculation, unit tests

include/numerics/type_vector.h

@@ -805,10 +805,11 @@ typename std::enable_if<
   TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
 TypeVector<T>::operator / (const Scalar & factor) const
 {
-  libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
-
   typedef typename CompareTypes<T, Scalar>::supertype TS;
 
+  libmesh_assert_not_equal_to (static_cast<TS>(factor),
+                               static_cast<TS>(0.));
+
 #if LIBMESH_DIM == 1
   return TypeVector<TS>(_coords[0]/factor);
 #endif
@@ -1069,6 +1070,33 @@ T solid_angle(const TypeVector<T> & v01,
 
 
 
+// Compute the (possibly 3D) circumcenter of three non-collinear points
+template <typename T>
+inline
+TypeVector<T> circumcenter(const TypeVector<T> & p0,
+                           const TypeVector<T> & p1,
+                           const TypeVector<T> & p2)
+{
+  // smallest subchords in Edge3 order, but this works for any order
+  const TypeVector<T> e02 = p2 - p0;
+  const TypeVector<T> e21 = p1 - p2;
+
+  const TypeVector<T> scaled_vec = e02.norm_sq()*e21 + e21.norm_sq()*e02;
+#if LIBMESH_DIM == 3
+  const TypeVector<T> numerator =
+    (e02.cross(e21)).cross(scaled_vec);
+#else
+  const T e02_cross_e21_z = e02(0)*e21(1)-e02(1)*e21(0);
+  const TypeVector<T> numerator {-e02_cross_e21_z*scaled_vec(1),
+                         e02_cross_e21_z*scaled_vec(0)};
+#endif
+  const TypeVector<T> e2C =
+    numerator*T(0.5)/cross_norm_sq(e02,e21);
+
+  return p2 + e2C;
+}
+
+
 
 /**
  * Compute |b x c|^2 without creating the extra temporary produced by

tests/numerics/type_vector_test.h

@@ -30,12 +30,14 @@
   CPPUNIT_TEST( testVectorAddAssign );          \
   CPPUNIT_TEST( testVectorSubAssign );          \
                                                 \
+  CPPUNIT_TEST( testIsZero );                   \
+  CPPUNIT_TEST( testSolidAngle );               \
+  CPPUNIT_TEST( testCircumcenter );             \
+                                                \
   CPPUNIT_TEST( testNormBase );                 \
   CPPUNIT_TEST( testNormSqBase );               \
   CPPUNIT_TEST( testValueBase );                \
   CPPUNIT_TEST( testZeroBase );                 \
-  CPPUNIT_TEST( testIsZero );                   \
-  CPPUNIT_TEST( testSolidAngle );               \
                                                 \
   CPPUNIT_TEST( testEqualityBase );             \
   CPPUNIT_TEST( testInEqualityBase );           \
@@ -390,6 +392,31 @@ class TypeVectorTestBase : public CppUnit::TestCase {
 #endif
   }
 
+  void testCircumcenter()
+  {
+    LOG_UNIT_TEST;
+
+    const DerivedClass origin(0), e_x(1);
+    DerivedClass twoone(2); twoone(1) = 1;
+
+    auto cc1 = circumcenter(origin, e_x, twoone);
+    LIBMESH_ASSERT_NUMBERS_EQUAL(cc1(0), 0.5, TOLERANCE*TOLERANCE);
+    LIBMESH_ASSERT_NUMBERS_EQUAL(cc1(1), 1.5, TOLERANCE*TOLERANCE);
+
+#if LIBMESH_DIM > 2
+    LIBMESH_ASSERT_NUMBERS_EQUAL(cc1(2), 0, TOLERANCE*TOLERANCE);
+    const DerivedClass twozeroone(2,0,1);
+    auto cc2 = circumcenter(origin, e_x, twozeroone);
+    LIBMESH_ASSERT_NUMBERS_EQUAL(cc2(0), 0.5, TOLERANCE*TOLERANCE);
+    LIBMESH_ASSERT_NUMBERS_EQUAL(cc2(1), 0, TOLERANCE*TOLERANCE);
+    LIBMESH_ASSERT_NUMBERS_EQUAL(cc2(2), 1.5, TOLERANCE*TOLERANCE);
+    auto cc3 = circumcenter(e_x, twoone, twozeroone);
+    LIBMESH_ASSERT_NUMBERS_EQUAL(cc3(0), Real(5)/3, TOLERANCE*TOLERANCE);
+    LIBMESH_ASSERT_NUMBERS_EQUAL(cc3(1), Real(1)/3, TOLERANCE*TOLERANCE);
+    LIBMESH_ASSERT_NUMBERS_EQUAL(cc3(2), Real(1)/3, TOLERANCE*TOLERANCE);
+#endif
+  }
+
   void testNormBase()
   {
     LOG_UNIT_TEST;