Merge pull request #1115 from LLNL/feature/kweiss/preshaped-geometry

Support replacement rules for preshaped geometry

RELEASE-NOTES.md

@@ -109,6 +109,9 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
 - Multimat: `MultiMat::makeOtherRelation()` now runs on the GPU with an appropriately-set allocator ID.
 - Multimat: `MultiMat::setCellMatRel(counts, indices)` now runs on the GPU, and accepts GPU-side data.
 - Renames `ArrayView::spacing()` to `ArrayView::minStride()`.
+- Klee: A shape's geometry no longer needs a `path` field when its `format` is "none"
+- Quest: Shapes without geometry can participate in replacement rules for sample-based shaping. Volume
+fractions for the associated materials must be supplied before shaping.
 
 ###  Fixed
 - Fixed issues with CUDA build in CMake versions 3.14.5 and above. Now require CMake 3.18+

src/axom/klee/Geometry.hpp

@@ -82,6 +82,9 @@ class Geometry
    */
   const std::string &getPath() const { return m_path; }
 
+  /// Predicate that returns true when the shape has an associated geometry
+  bool hasGeometry() const { return !m_path.empty(); }
+
   /**
    * Get a GeometryOperator to apply to this geometry. Can be null.
    *

src/axom/klee/IO.cpp

@@ -73,8 +73,8 @@ struct FromInlet<axom::klee::GeometryData>
   axom::klee::GeometryData operator()(const axom::inlet::Container &base)
   {
     axom::klee::GeometryData data;
-    data.format = base["format"];
-    data.path = base["path"];
+    data.format = base.contains("format") ? base.get<std::string>("format") : "";
+    data.path = base.contains("path") ? base.get<std::string>("path") : "";
     data.operatorData =
       base["operators"].get<axom::klee::internal::GeometryOperatorData>();
 
@@ -116,9 +116,9 @@ using inlet::Inlet;
 void defineGeometry(Container &geometry)
 {
   geometry.addString("format", "The format of the input file").required();
-  geometry
-    .addString("path", "The path of the input file, relative to the yaml file")
-    .required();
+  geometry.addString("path",
+                     "The path of the input file, relative to the yaml file."
+                     "Required unless 'format' is 'none'");
   internal::defineDimensionsField(
     geometry,
     "start_dimensions",
@@ -150,6 +150,11 @@ void defineShapeList(Inlet &document)
                            "The list of materials this shape replaces");
   shapeList.addStringArray("does_not_replace",
                            "The list of materials this shape does not replace");
+  auto &geometry =
+    shapeList.addStruct("geometry",
+                        "Contains information about the shape's geometry");
+  defineGeometry(geometry);
+
   // Verify syntax here, semantics later!!!
   shapeList.registerVerifier(
     [](const Container &shape,
@@ -162,12 +167,25 @@ void defineShapeList(Inlet &document)
           errors);
         return false;
       }
+
+      if(shape.contains("geometry"))
+      {
+        const auto geom = shape.get<GeometryData>("geometry");
+        if(geom.path.empty() && geom.format != "none")
+        {
+          INLET_VERIFICATION_WARNING(
+            shape.name(),
+            fmt::format("'geometry/path' field required unless "
+                        "'geometry/format' is 'none'. "
+                        "Provided format was '{}'",
+                        geom.format),
+            errors);
+          return false;
+        }
+      }
+
       return true;
     });
-  auto &geometry =
-    shapeList.addStruct("geometry",
-                        "Contains information about the shape's geometry");
-  defineGeometry(geometry);
 }
 
 /**

src/axom/klee/tests/klee_geometry.cpp

@@ -25,6 +25,8 @@ TEST(GeometryTest, dimensions_noOperators)
   Geometry geometry {startProperties, "test format", "test path", nullptr};
   EXPECT_EQ(startProperties, geometry.getStartProperties());
   EXPECT_EQ(startProperties, geometry.getEndProperties());
+
+  EXPECT_TRUE(geometry.hasGeometry());
 }
 
 TEST(GeometryTest, dimensions_dimensionPreservingOperator)
@@ -43,5 +45,14 @@ TEST(GeometryTest, dimensions_dimensionPreservingOperator)
   EXPECT_EQ(endProperties, geometry.getEndProperties());
 }
 
+TEST(GeometryTest, emptyPath)
+{
+  TransformableGeometryProperties startProperties {Dimensions::Three,
+                                                   LengthUnit::mils};
+  Geometry geometry {startProperties, "none", "", nullptr};
+
+  EXPECT_FALSE(geometry.hasGeometry());
+}
+
 }  // namespace klee
 }  // namespace axom

src/axom/klee/tests/klee_io.cpp

@@ -173,24 +173,51 @@ TEST(IOTest, readShapeSet_missingMaterial)
 
 TEST(IOTest, readShapeSet_missingGeometryPath)
 {
-  auto input = R"(
-    dimensions: 2
-
-    shapes:
-      - name: wheel
-        material: steel
-        geometry:
-          format: test_format
-  )";
-
-  try
+  // Expect a validation error when 'geometry/path' is missing
+  // and 'geometry/format' is not "none"
   {
-    readShapeSetFromString(input);
-    FAIL() << "Should have thrown";
+    auto input = R"(
+      dimensions: 2
+
+      shapes:
+        - name: wheel
+          material: steel
+          geometry:
+            format: test_format
+    )";
+
+    try
+    {
+      readShapeSetFromString(input);
+      FAIL() << "Should have thrown";
+    }
+    catch(const KleeError &err)
+    {
+      EXPECT_THAT(err.what(), HasSubstr("Provided format"));
+    }
   }
-  catch(const KleeError &err)
+
+  // Valid when 'geometry/path' is missing and 'geometry/format' is none
   {
-    EXPECT_THAT(err.what(), HasSubstr("path"));
+    auto input = R"(
+      dimensions: 2
+
+      shapes:
+        - name: wheel
+          material: steel
+          geometry:
+            format: none
+    )";
+
+    try
+    {
+      readShapeSetFromString(input);
+      SUCCEED();
+    }
+    catch(const KleeError &err)
+    {
+      FAIL() << "Should not have thrown. Error message: " << err.what();
+    }
   }
 }
 

src/axom/multimat/multimat.cpp

@@ -488,7 +488,6 @@ void ScanRelationOffsetsRAJA(const axom::ArrayView<const IndexType> counts,
                              const axom::ArrayView<IndexType> firstIndices)
 {
 #if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE)
-  using AtomicPolicy = typename axom::execution_space<ExecSpace>::atomic_policy;
   using LoopPolicy = typename axom::execution_space<ExecSpace>::loop_policy;
 
   // The begins array has one more element than the counts array. The last

src/axom/primal/operators/winding_number.hpp

@@ -268,9 +268,13 @@ double winding_number(const Point<T, 3>& q,
                       const double edge_tol = 1e-8,
                       const double EPS = 1e-8)
 {
+  using Vec3 = Vector<T, 3>;
+
   if(tri.area() == 0) return 0;
 
-  Vector<T, 3> a(q, tri[0]), b(q, tri[1]), c(q, tri[2]);
+  const Vec3 a = tri[0] - q;
+  const Vec3 b = tri[1] - q;
+  const Vec3 c = tri[2] - q;
 
   // Compute norms. Possibly return early
   const double a_norm = a.norm();
@@ -279,40 +283,34 @@ double winding_number(const Point<T, 3>& q,
 
   if(a_norm < edge_tol || b_norm < edge_tol || c_norm < edge_tol) return 0;
 
-  double num = Vector<T, 3>::scalar_triple_product(a, b, c);
+  const double num = Vec3::scalar_triple_product(a, b, c);
   if(axom::utilities::isNearlyEqual(num, 0.0, EPS))
   {
     isOnFace = true;
     return 0;
   }
 
-  double denom = a_norm * b_norm * c_norm +
-    a_norm * Vector<T, 3>::dot_product(b, c) +
-    b_norm * Vector<T, 3>::dot_product(a, c) +
-    c_norm * Vector<T, 3>::dot_product(a, b);
+  const double denom = a_norm * b_norm * c_norm  //
+    + a_norm * b.dot(c) + b_norm * a.dot(c) + c_norm * a.dot(b);
 
   // Handle direct cases where argument to atan is undefined
   if(axom::utilities::isNearlyEqual(denom, 0.0, EPS))
   {
     return (num > 0) ? 0.25 : -0.25;
   }
 
+  // Note: denom==0 and num==0 handled above
   if(denom > 0)
   {
     return 0.5 * M_1_PI * atan(num / denom);
   }
   else
   {
-    if(num > 0)
-    {
-      return 0.5 * M_1_PI * atan(num / denom) + 0.5;
-    }
-    if(num < 0)
-    {
-      return 0.5 * M_1_PI * atan(num / denom) - 0.5;
-    }
+    return (num > 0) ? 0.5 * M_1_PI * atan(num / denom) + 0.5
+                     : 0.5 * M_1_PI * atan(num / denom) - 0.5;
   }
 }
+
 /*!
  * \brief Computes the solid angle winding number for a 3D triangle
  *

src/axom/primal/tests/primal_bezier_patch.cpp

@@ -824,7 +824,6 @@ TEST(primal_bezierpatch, rational_evaluation_split)
   const int DIM = 3;
   using CoordType = double;
   using PointType = primal::Point<CoordType, DIM>;
-  using VectorType = primal::Vector<CoordType, DIM>;
   using BezierPatchType = primal::BezierPatch<CoordType>;
 
   const int ord_u = 3;

src/axom/primal/tests/primal_solid_angle.cpp

@@ -22,7 +22,6 @@ namespace primal = axom::primal;
 TEST(primal_solid_angle, triangle)
 {
   using Point3D = primal::Point<double, 3>;
-  using Vector3D = primal::Vector<double, 3>;
   using Triangle = primal::Triangle<double, 3>;
 
   Point3D origin {0.0, 0.0, 0.0};
@@ -207,7 +206,6 @@ TEST(primal_solid_angle, degenerate_polygon)
 {
   using Point3D = primal::Point<double, 3>;
   using Vector3D = primal::Vector<double, 3>;
-  using Triangle = primal::Triangle<double, 3>;
   using Polygon = primal::Polygon<double, 3>;
 
   Vector3D v1 = Vector3D({0.0, 1.0, 2.0}).unitVector();
@@ -282,7 +280,6 @@ TEST(primal_solid_angle, selfintersecting_star)
 {
   using Point3D = primal::Point<double, 3>;
   using Vector3D = primal::Vector<double, 3>;
-  using Triangle = primal::Triangle<double, 3>;
   using Polygon = primal::Polygon<double, 3>;
 
   Point3D single_query {1, 2, 3};
@@ -383,7 +380,6 @@ TEST(primal_solid_angle, selfintersecting_star)
 TEST(primal_solid_angle, selfintersecting_quadrilateral)
 {
   using Point3D = primal::Point<double, 3>;
-  using Vector3D = primal::Vector<double, 3>;
   using Triangle = primal::Triangle<double, 3>;
   using Polygon = primal::Polygon<double, 3>;