Feature/mir

.gitignore

@@ -14,3 +14,4 @@ _axom_build_and_test_*
 tpl_dirs_summary.json
 *.swp
 *.vscode*
+uberenv_libs

src/axom/CMakeLists.txt

@@ -30,6 +30,7 @@ axom_add_component(COMPONENT_NAME primal DEFAULT_STATE ${AXOM_ENABLE_ALL_COMPONE
 axom_add_component(COMPONENT_NAME spin   DEFAULT_STATE ${AXOM_ENABLE_ALL_COMPONENTS})
 axom_add_component(COMPONENT_NAME sidre  DEFAULT_STATE ${AXOM_ENABLE_ALL_COMPONENTS})
 axom_add_component(COMPONENT_NAME quest  DEFAULT_STATE ${AXOM_ENABLE_ALL_COMPONENTS})
+axom_add_component(COMPONENT_NAME mir    DEFAULT_STATE ${AXOM_ENABLE_ALL_COMPONENTS})
 
 # Combine all component object libraries into a unified library
 blt_add_library(NAME       axom

src/axom/config.hpp.in

@@ -90,6 +90,7 @@
  * Compiler defines for the toolkit components
  */
 #cmakedefine AXOM_USE_MINT
+#cmakedefine AXOM_USE_MIR
 #cmakedefine AXOM_USE_LUMBERJACK
 #cmakedefine AXOM_USE_PRIMAL
 #cmakedefine AXOM_USE_QUEST

src/axom/core/tests/utils_utilities.cpp

@@ -116,3 +116,34 @@ TEST(core_Utilities,minmax)
     EXPECT_EQ(a, axom::utilities::max(a,b));
   }
 }
+
+TEST(core_Utilities, lerp)
+{
+  std::cout<<"Testing linear interpolation (lerp) function."<< std::endl;
+
+  double f0 = 50.0;
+  double f1 = 100.0;
+
+  // Test end points
+  {
+    EXPECT_DOUBLE_EQ( f0, axom::utilities::lerp(f0, f1, 0.) );
+    EXPECT_DOUBLE_EQ( f1, axom::utilities::lerp(f1, f0, 0.) );
+
+    EXPECT_DOUBLE_EQ( f1, axom::utilities::lerp(f0, f1, 1.) );
+    EXPECT_DOUBLE_EQ( f0, axom::utilities::lerp(f1, f0, 1.) );
+  }
+
+  // Test midpoint
+  {
+    double t = 0.5;
+    double exp = 75.;
+    EXPECT_DOUBLE_EQ(  exp, axom::utilities::lerp(f0, f1, t));
+  }
+
+  // Another test
+  {
+    double t = 0.66;
+    double exp = 83.;
+    EXPECT_DOUBLE_EQ(  exp, axom::utilities::lerp(f0, f1, t));
+  }
+}

src/axom/core/utilities/Utilities.hpp

@@ -96,6 +96,22 @@ inline T log2( T& val)
   return std::log2(val);
 }
 
+
+/*!
+ * \brief Linearly interpolates between two values
+ * \param [in] val0 The first value
+ * \param [in] val2 The second value
+ * \param [in] t The interpolation parameter.
+ * \return The interpolated value
+ */
+template < typename T >
+inline AXOM_HOST_DEVICE
+T lerp( T v0, T v1, T t)
+{
+  constexpr T one = T(1);
+  return (one-t)*v0 + t*v1;
+}
+
 /*!
  * \brief Clamps an input value to a given range.
  * \param [in] val  The value to clamp.
@@ -113,7 +129,6 @@ T clampVal( T val, T lower, T upper )
          : (val > upper) ? upper : val;
 }
 
-
 /*!
  * \brief Clamps the upper range on an input value
  *

src/axom/mainpage.md

@@ -8,6 +8,7 @@ Axom provides libraries that address common computer science needs.  It grew fro
 * @subpage coretop provides shared utility functionality to all components.
 * @subpage lumberjacktop provides logging aggregation and filtering capability.
 * @subpage minttop provides a comprehensive mesh data model.
+* @subpage mirtop provides algorithms for material interface reconstruction on multimaterial meshes.
 * @subpage primaltop provides an API for geometric primitives and computational geometry tests.
 * @subpage questtop provides an API to query point distance and position relative to meshes.
 * @subpage sidretop provides a data store with hierarchical structure.
@@ -20,6 +21,7 @@ Dependencies between components are as follows:
 - Slic optionally depends on Lumberjack
 - Slam, Primal, Mint, Quest, Spin, and Sidre depend on Slic
 - Mint optionally depends on Sidre
+- Mir depends on Slic, Slam and Primal
 - Spin depends on Primal and Slam
 - Quest depends on Slam, Primal, Spin, and Mint
 

src/axom/mir/CMakeLists.txt

@@ -0,0 +1,76 @@
+# Copyright (c) 2017-2019, Lawrence Livermore National Security, LLC and
+# other Axom Project Developers. See the top-level COPYRIGHT file for details.
+#
+# SPDX-License-Identifier: (BSD-3-Clause)
+#------------------------------------------------------------------------------
+# MIR -- Material Interface Reconstruction
+#------------------------------------------------------------------------------
+
+#------------------------------------------------------------------------------
+# Check necessary dependencies
+#------------------------------------------------------------------------------
+axom_component_requires(NAME       MIR
+                        COMPONENTS SLIC SLAM PRIMAL)
+
+#------------------------------------------------------------------------------
+# Specify all headers/sources
+#------------------------------------------------------------------------------
+set(mir_headers
+    MIRMesh.hpp
+    MIRMeshTypes.hpp
+    ZooClippingTables.hpp
+    InterfaceReconstructor.hpp
+    CellData.hpp
+    MeshTester.hpp
+    CellClipper.hpp
+    MIRUtilities.hpp
+    CellGenerator.hpp
+    )
+
+set(mir_sources
+    MIRMesh.cpp
+    InterfaceReconstructor.cpp
+    ZooClippingTables.cpp
+    CellData.cpp
+    MeshTester.cpp
+    CellClipper.cpp
+    CellGenerator.cpp
+    )
+
+#------------------------------------------------------------------------------
+# Build and install the library
+#------------------------------------------------------------------------------
+set(mir_depends_on core slic slam)
+
+blt_add_library(
+    NAME        mir
+    SOURCES     ${mir_sources}
+    HEADERS     ${mir_headers}
+    DEPENDS_ON  ${mir_depends_on}
+    FOLDER      axom/mir
+    OBJECT      TRUE )
+
+axom_write_unified_header(NAME    mir
+                          HEADERS ${mir_headers} )
+
+axom_install_component(NAME    mir
+                       HEADERS ${mir_headers} )
+
+#------------------------------------------------------------------------------
+# Add tests, benchmarks and examples
+#------------------------------------------------------------------------------
+if (AXOM_ENABLE_TESTS)
+  add_subdirectory(tests)
+endif()
+
+if (AXOM_ENABLE_EXAMPLES)
+  add_subdirectory(examples)
+endif()
+
+
+#------------------------------------------------------------------------------
+# Add code checks
+#------------------------------------------------------------------------------
+axom_add_code_checks(
+        PREFIX mir )
+

src/axom/mir/CellClipper.cpp

@@ -0,0 +1,157 @@
+// Copyright (c) 2017-2019, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level COPYRIGHT file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#include "CellClipper.hpp"
+
+namespace axom
+{
+namespace mir
+{
+
+//--------------------------------------------------------------------------------
+
+CellClipper::CellClipper()
+{
+
+}
+
+//--------------------------------------------------------------------------------
+
+CellClipper::~CellClipper()
+{
+
+}
+
+//--------------------------------------------------------------------------------
+
+// Computes the t-values where each edge is clipped, as well as the topology of the new output cells after clipping the original cell
+// Outputs the newElements, newVertices maps and the verticesClippingTValue array[]
+void CellClipper::computeClippingPoints(const mir::Shape shapeType, 
+                                        const std::vector<std::vector<axom::float64> >& vertexVF, 
+                                        std::map<int, std::vector<int> >& newElements, 
+                                        std::map<int, std::vector<int> >& newVertices, 
+                                        axom::float64* tValues)
+{
+  // Determine the clipping case for the current element
+  unsigned int caseIndex = determineClippingCase( shapeType, vertexVF[0], vertexVF[1] );
+
+  std::vector<std::vector<int> > clipTable = getClipTable(shapeType);
+
+  // Create the new polygons based on the clipping case
+  int currentElementIndex = 0;    // the next available element index
+  int i = 0;
+  int numVertices = clipTable[caseIndex][i];
+
+  // for each new element in the current clipping case
+  while (numVertices != -1)   
+  {
+    // for each vertex of the new element
+    for (int j = 0; j < numVertices; ++j)
+    {
+      // Find the id of the next vertex of the new element
+      int vID = clipTable[caseIndex][i + (j+1)];
+
+      // Associate the vertex and element together
+      newElements[currentElementIndex].push_back(vID);
+      newVertices[vID].push_back(currentElementIndex);
+
+      if ( vID >= mir::utilities::numVerts(shapeType) )
+      {
+        int vertexOneID = mir::utilities::getEdgeEndpoint(shapeType, vID, true);
+        int vertexTwoID = mir::utilities::getEdgeEndpoint(shapeType, vID, false);
+
+        tValues[vID] = computeTValueOnEdge( vertexVF[0][vertexOneID], vertexVF[1][vertexOneID], vertexVF[0][vertexTwoID], vertexVF[1][vertexTwoID] );
+      }
+    }
+
+    // Increment the element index counter, marking the current element as being finished processed
+    currentElementIndex++;
+
+    // Increase index into lookup table to the next element
+    i += (numVertices + 1);
+    numVertices = clipTable[caseIndex][i];
+  }
+}
+
+//--------------------------------------------------------------------------------
+
+unsigned int CellClipper::determineClippingCase(const mir::Shape shapeType, 
+                                                const std::vector<axom::float64>& matOneVF, 
+                                                const std::vector<axom::float64>& matTwoVF)
+{
+  unsigned int caseIndex = 0;
+
+  int numVertices = mir::utilities::numVerts(shapeType);
+
+  for (int vID = 0; vID < numVertices; ++vID)
+  {
+    if (matOneVF[vID] > matTwoVF[vID])
+    {
+      unsigned int bitIndex = (numVertices - 1) - vID;
+
+      caseIndex |= (1 << bitIndex);
+    } 
+  }
+
+  return caseIndex;
+}
+
+//--------------------------------------------------------------------------------
+
+axom::float64  CellClipper::computeTValueOnEdge(axom::float64 vfMatOneVertexOne,
+                                                axom::float64 vfMatTwoVertexOne, 
+                                                axom::float64 vfMatOneVertexTwo, 
+                                                axom::float64 vfMatTwoVertexTwo)
+{
+  axom::float64 ret = 0.0;
+
+  // TODO: Perhaps just handle NULL_MAT by return 0, since that is what will happen anyways?
+  // Handle NULL_MAT, which has a vf of -1.0, but which needs to be 0.0 for the purposes of computing the clipping point
+  if (vfMatOneVertexOne < 0.0) { vfMatOneVertexOne = 0.0; };
+  if (vfMatTwoVertexOne < 0.0) { vfMatTwoVertexOne = 0.0; };
+  if (vfMatOneVertexTwo < 0.0) { vfMatOneVertexTwo = 0.0; };
+  if (vfMatTwoVertexTwo < 0.0) { vfMatTwoVertexTwo = 0.0; };
+
+  axom::float64 numerator = vfMatTwoVertexOne - vfMatOneVertexOne;
+  axom::float64 denominator = -vfMatOneVertexOne + vfMatOneVertexTwo + vfMatTwoVertexOne - vfMatTwoVertexTwo;
+
+  if (denominator != 0.0)
+  {
+    ret = numerator / denominator;
+  }
+
+  if (ret > 1.0 || ret < 0.0)
+  {
+    // This shouldn't happen...
+    printf("    OUT OF BOUNDS T VALUE: %f\n", ret);
+
+    // Clamp the t value
+    ret = fmin(1.0, ret);
+    ret = fmax(0.0, ret);
+  }
+
+  return ret;
+}
+
+//--------------------------------------------------------------------------------
+
+const std::vector<std::vector<int> >& CellClipper::getClipTable(const mir::Shape shapeType)
+{
+  switch ( shapeType )
+  {
+    case mir::Shape::Triangle:
+      return triangleClipTableVec;
+    case mir::Shape::Quad:
+      return quadClipTableVec;
+    default:
+      printf("No clipping table for this shape type.\n");
+      return triangleClipTableVec;
+  }
+}
+
+//--------------------------------------------------------------------------------
+
+}
+}