From 4b44dd52b25b797710264c48f36bafad009c7336 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 19 May 2026 10:40:39 -0500 Subject: [PATCH 1/4] Add circumcenter calculation, unit tests --- include/numerics/type_vector.h | 27 +++++++++++++++++++++++++++ tests/numerics/type_vector_test.h | 30 ++++++++++++++++++++++++++++-- 2 files changed, 55 insertions(+), 2 deletions(-) diff --git a/include/numerics/type_vector.h b/include/numerics/type_vector.h index aaf79a9fd22..783631cd0a7 100644 --- a/include/numerics/type_vector.h +++ b/include/numerics/type_vector.h @@ -1069,6 +1069,33 @@ T solid_angle(const TypeVector & v01, +// Compute the (possibly 3D) circumcenter of three non-collinear points +template +inline +TypeVector circumcenter(const TypeVector & p0, + const TypeVector & p1, + const TypeVector & p2) +{ + // smallest subchords in Edge3 order, but this works for any order + const TypeVector e02 = p2 - p0; + const TypeVector e21 = p1 - p2; + + const TypeVector scaled_vec = e02.norm_sq()*e21 + e21.norm_sq()*e02; +#if LIBMESH_DIM == 3 + const TypeVector numerator = + (e02.cross(e21)).cross(scaled_vec); +#else + const T e02_cross_e21_z = e02(0)*e21(1)-e02(1)*e21(0); + const TypeVector numerator {-e02_cross_e21_z*scaled_vec(1), + e02_cross_e21_z*scaled_vec(0)}; +#endif + const TypeVector e2C = + numerator/2/cross_norm_sq(e02,e21); + + return p2 + e2C; +} + + /** * Compute |b x c|^2 without creating the extra temporary produced by diff --git a/tests/numerics/type_vector_test.h b/tests/numerics/type_vector_test.h index ede8abc37b0..12126246646 100644 --- a/tests/numerics/type_vector_test.h +++ b/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,30 @@ class TypeVectorTestBase : public CppUnit::TestCase { #endif } + void testCircumcenter() + { + LOG_UNIT_TEST; + + const DerivedClass origin(0), e_x(1), twoone(2,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; From 9a6c7f891a4efb30fe41409fdcc69e86b5d76c16 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 19 May 2026 14:04:45 -0500 Subject: [PATCH 2/4] Do assert with supertype The assert itself should have been handling conversions, but for complex/int the direct conversion wasn't allowed. --- include/numerics/type_vector.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/numerics/type_vector.h b/include/numerics/type_vector.h index 783631cd0a7..f340d6a6061 100644 --- a/include/numerics/type_vector.h +++ b/include/numerics/type_vector.h @@ -805,10 +805,11 @@ typename std::enable_if< TypeVector::supertype>>::type TypeVector::operator / (const Scalar & factor) const { - libmesh_assert_not_equal_to (factor, static_cast(0.)); - typedef typename CompareTypes::supertype TS; + libmesh_assert_not_equal_to (static_cast(factor), + static_cast(0.)); + #if LIBMESH_DIM == 1 return TypeVector(_coords[0]/factor); #endif From 5a6689cc8bed3a3912369bbe57a050ffbe503669 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 19 May 2026 14:06:52 -0500 Subject: [PATCH 3/4] Multiply by T(0.5), don't divide by 2 This fixes T=complex issues and is more efficient if the compiler doesn't inline here. --- include/numerics/type_vector.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/numerics/type_vector.h b/include/numerics/type_vector.h index f340d6a6061..00efee3fa61 100644 --- a/include/numerics/type_vector.h +++ b/include/numerics/type_vector.h @@ -1091,7 +1091,7 @@ TypeVector circumcenter(const TypeVector & p0, e02_cross_e21_z*scaled_vec(0)}; #endif const TypeVector e2C = - numerator/2/cross_norm_sq(e02,e21); + numerator*T(0.5)/cross_norm_sq(e02,e21); return p2 + e2C; } From ffc20ed09ee00e10aaf875df87331c07fafcf032 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 19 May 2026 14:07:36 -0500 Subject: [PATCH 4/4] Initialize vec --enable-complex compatibly We really need to fix this, not keep working around it... --- tests/numerics/type_vector_test.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/numerics/type_vector_test.h b/tests/numerics/type_vector_test.h index 12126246646..1a8e6afb8ad 100644 --- a/tests/numerics/type_vector_test.h +++ b/tests/numerics/type_vector_test.h @@ -396,7 +396,8 @@ class TypeVectorTestBase : public CppUnit::TestCase { { LOG_UNIT_TEST; - const DerivedClass origin(0), e_x(1), twoone(2,1); + 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);