diff --git a/.gitignore b/.gitignore index 501c8041a6..5ab2039dfa 100644 --- a/.gitignore +++ b/.gitignore @@ -14,3 +14,4 @@ _axom_build_and_test_* tpl_dirs_summary.json *.swp *.vscode* +uberenv_libs diff --git a/src/axom/CMakeLists.txt b/src/axom/CMakeLists.txt index 4dc4d676d4..467e1d0f7c 100644 --- a/src/axom/CMakeLists.txt +++ b/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 diff --git a/src/axom/config.hpp.in b/src/axom/config.hpp.in index 9992a767c1..399d296c6b 100644 --- a/src/axom/config.hpp.in +++ b/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 diff --git a/src/axom/core/tests/utils_utilities.cpp b/src/axom/core/tests/utils_utilities.cpp index 85c092b7e8..1b11d57a25 100644 --- a/src/axom/core/tests/utils_utilities.cpp +++ b/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)); + } +} diff --git a/src/axom/core/utilities/Utilities.hpp b/src/axom/core/utilities/Utilities.hpp index d891db13e2..522cea8509 100644 --- a/src/axom/core/utilities/Utilities.hpp +++ b/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 * diff --git a/src/axom/mainpage.md b/src/axom/mainpage.md index 827988937a..366e41642b 100644 --- a/src/axom/mainpage.md +++ b/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 diff --git a/src/axom/mir/CMakeLists.txt b/src/axom/mir/CMakeLists.txt new file mode 100644 index 0000000000..ef829028f0 --- /dev/null +++ b/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 ) + diff --git a/src/axom/mir/CellClipper.cpp b/src/axom/mir/CellClipper.cpp new file mode 100644 index 0000000000..c5f25d4272 --- /dev/null +++ b/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 >& vertexVF, + std::map >& newElements, + std::map >& newVertices, + axom::float64* tValues) +{ + // Determine the clipping case for the current element + unsigned int caseIndex = determineClippingCase( shapeType, vertexVF[0], vertexVF[1] ); + + std::vector > 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& matOneVF, + const std::vector& 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 >& 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; + } +} + +//-------------------------------------------------------------------------------- + +} +} \ No newline at end of file diff --git a/src/axom/mir/CellClipper.hpp b/src/axom/mir/CellClipper.hpp new file mode 100644 index 0000000000..fcdd650db3 --- /dev/null +++ b/src/axom/mir/CellClipper.hpp @@ -0,0 +1,120 @@ +// 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) + +/** + * \file CellClipper.hpp + * + * \brief Contains the specification for the CellClipper class. + * + */ + +#ifndef __CELL_CLIPPER_H +#define __CELL_CLIPPER_H + +#include "axom/core.hpp" // for axom macros +#include "axom/slam.hpp" // unified header for slam classes and functions + +#include "MIRMesh.hpp" +#include "MIRUtilities.hpp" +#include "MIRMeshTypes.hpp" +#include "CellData.hpp" +#include "ZooClippingTables.hpp" +#include "MIRUtilities.hpp" + +//-------------------------------------------------------------------------------- + +namespace axom +{ +namespace mir +{ + +//-------------------------------------------------------------------------------- + + /** + * \class CellClipper + * + * \brief A class that contains the functionality for taking an input cell + * and determining how it should be clipped. + * + */ + class CellClipper + { + public: + /** + * \brief Default constructor. + */ + CellClipper(); + + /** + * \brief Default destructor. + */ + ~CellClipper(); + + /** + * \brief Computes the elements, vertices, and t values resulting from splitting the element with the given vertex volume fractions. + * + * \param shapeType The shape type of the element. + * \param vertexVF The vertex volume fractions of the two materials with which to clip the cell. + * \param newElements An ordered map of the generated elements' IDs to a list of the vertex IDs in the local frame of the output element. + * \param newVertices An order map of the vertex IDs in the local frame of the output element to the generated elements' IDs it is associated with. + * \param tValues An array of t values where each of the midpoint vertices are for the newly generated elements. + * + */ + void computeClippingPoints(const mir::Shape shapeType, + const std::vector >& vertexVF, + std::map >& newElements, + std::map >& newVertices, + axom::float64* tValues); + + /** + * \brief Determines the index into the clipping table of the given shape based on the volume fractions at its vertices. + * + * \param shapeType The shape type of the element. + * \param matOneVF The volume fractions at the vertices of the element for the first material. + * \param matTwoVF The volume fractions at the vertices of the element for the second material. + * + * \return The index into the clipping table. + */ + unsigned int determineClippingCase(const mir::Shape shapeType, + const std::vector& matOneVF, + const std::vector& matTwoVF); + + + /** + * \brief Computes the t value as a percent from vertex one to vertex two based on the materials given. + * + * \param vfMatOneVertexOne The volume fraction of material one present at vertex one. + * \param vfMatTwoVertexOne The volume fraction of material two present at vertex one. + * \param vfMatOneVertexTwo The volume fraction of material one present at vertex two. + * \param vfMatTwoVertexTwo The volume fraction of material two present at vertex two. + * + * \return The percent of the distance from vertex one to vertex two where the edge should be clipped. + * + * \note When one material's volume fractions dominates the other material's, then the edge should not be clipped and this function will return 0.0. + * \note When one of the materials is the NULL_MAT (and has vf = -1.0), these values are set to 0 in order to interpolate properly. + */ + axom::float64 computeTValueOnEdge(axom::float64 vfMatOneVertexOne, + axom::float64 vfMatTwoVertexOne, + axom::float64 vfMatOneVertexTwo, + axom::float64 vfMatTwoVertexTwo); + + private: + /** + * \brief Returns a reference to the appropriate clipping table to use for the shape type. + * + * \param The shape type of the element. + * + * \return A reference to the clipping table. + */ + const std::vector >& getClipTable(const mir::Shape shapeType); + + }; + +//-------------------------------------------------------------------------------- + +} +} + +#endif \ No newline at end of file diff --git a/src/axom/mir/CellData.cpp b/src/axom/mir/CellData.cpp new file mode 100644 index 0000000000..176d1d6fc7 --- /dev/null +++ b/src/axom/mir/CellData.cpp @@ -0,0 +1,100 @@ +// 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 "CellData.hpp" + +namespace axom +{ +namespace mir +{ + + //-------------------------------------------------------------------------------- + + CellData::CellData() + { + + } + + //-------------------------------------------------------------------------------- + + CellData::~CellData() + { + + } + + //-------------------------------------------------------------------------------- + + void CellData::mergeCell(const CellData& cellToMerge) + { + // Initialize index offsets + int evBeginsOffset = m_topology.m_evInds.size(); + int veBeginsOffset = m_topology.m_veInds.size(); + + int vertexIndexOffset = m_numVerts; + int elementIndexOffset = m_numElems; + + // Merge the cell topology information + for (unsigned long i = 0; i < cellToMerge.m_topology.m_evInds.size(); ++i) + { + m_topology.m_evInds.push_back(cellToMerge.m_topology.m_evInds[i] + vertexIndexOffset); + } + + for (unsigned long i = 1; i < cellToMerge.m_topology.m_evBegins.size(); ++i) + { + m_topology.m_evBegins.push_back(cellToMerge.m_topology.m_evBegins[i] + evBeginsOffset); + } + + for (unsigned long i = 0; i < cellToMerge.m_topology.m_veInds.size(); ++i) + { + m_topology.m_veInds.push_back(cellToMerge.m_topology.m_veInds[i] + elementIndexOffset); + } + + for (unsigned long i = 1; i < cellToMerge.m_topology.m_veBegins.size(); ++i) + { + m_topology.m_veBegins.push_back(cellToMerge.m_topology.m_veBegins[i] + veBeginsOffset); + } + + // Merge the vertex positions + for (unsigned long i = 0; i < cellToMerge.m_mapData.m_vertexPositions.size(); ++i) + { + m_mapData.m_vertexPositions.push_back(cellToMerge.m_mapData.m_vertexPositions[i]); + } + + // Merge the vertex volume fractions + for (unsigned long matID = 0; matID < m_mapData.m_vertexVolumeFractions.size(); ++matID) + { + for (unsigned long vID = 0; vID < cellToMerge.m_mapData.m_vertexVolumeFractions[matID].size(); ++vID) + { + m_mapData.m_vertexVolumeFractions[matID].push_back(cellToMerge.m_mapData.m_vertexVolumeFractions[matID][vID]); + } + } + + // Merge the elements' dominant materials + for (unsigned long i = 0; i < cellToMerge.m_mapData.m_elementDominantMaterials.size(); ++i) + { + m_mapData.m_elementDominantMaterials.push_back(cellToMerge.m_mapData.m_elementDominantMaterials[i]); + } + + // Merge the elements' parent ids + for (unsigned long i = 0; i < cellToMerge.m_mapData.m_elementParents.size(); ++i) + { + m_mapData.m_elementParents.push_back(cellToMerge.m_mapData.m_elementParents[i]); + } + + // Merge the elements' shape types + for (unsigned long i = 0; i < cellToMerge.m_mapData.m_shapeTypes.size(); ++i) + { + m_mapData.m_shapeTypes.push_back(cellToMerge.m_mapData.m_shapeTypes[i]); + } + + // Merge the total number of verts and elems in the resulting cell + m_numVerts += cellToMerge.m_numVerts; + m_numElems += cellToMerge.m_numElems; + } + + //-------------------------------------------------------------------------------- + +} +} \ No newline at end of file diff --git a/src/axom/mir/CellData.hpp b/src/axom/mir/CellData.hpp new file mode 100644 index 0000000000..a73134374c --- /dev/null +++ b/src/axom/mir/CellData.hpp @@ -0,0 +1,97 @@ +// 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) + +/** + * \file CellData.hpp + * + * \brief Contains the specifications for the CellData class + * and CellTopologyData and CellMapData structs. + */ + +#ifndef __CELL_DATA_H__ +#define __CELL_DATA_H__ + +#include "axom/core.hpp" // for axom macros +#include "axom/slam.hpp" + +#include "MIRMeshTypes.hpp" + +namespace numerics = axom::numerics; +namespace slam = axom::slam; + +namespace axom +{ +namespace mir +{ + + /** + * \struct CellTopologyData + * + * \brief Struct for collecting data that specifies a mesh or cell's connectivity/topology. + */ + struct CellTopologyData + { + std::vector m_evInds; + std::vector m_evBegins; + std::vector m_veInds; + std::vector m_veBegins; + }; + + /** + * \struct CellMapData + * + * \brief Struct for collecting the data that will populate a mesh's map data structures. + */ + struct CellMapData + { + std::vector m_vertexPositions; // Data that goes into MIRMesh's vertexPositions PointMap + std::vector > m_vertexVolumeFractions; // Data that goes into MIRMesh's materialVolumeFractionsVertex ScalarMap + std::vector m_elementDominantMaterials; // Data that goes into MIRMesh's elementDominantColors IntMap + std::vector m_elementParents; // Data that goes into MIRMesh's elementParentIDs IntMap + std::vector m_shapeTypes; // Data that goes into MIRMesh's shapeType IntMap + }; + + /** + * \class CellData + * + * \brief The CellData class represents an arbitrary number of cells that are + * within the same local coordinate system (i.e. share a set of vertices + * and elements). + * + * \detail This class is intended to be used as a helper class to hold intermediate + * data while processing a mesh, and to be used as input to the MIRMesh class + * to fully initialize it. + */ + class CellData + { + public: + /** + * \brief Default constructor. + */ + CellData(); + + /** + * \brief Default destructor. + */ + ~CellData(); + + /** + * \brief Merges the cell data from the given cell into this cell. + * + * \param cellToMerge The cell whose data will be taken and merged in. + */ + void mergeCell(const CellData& cellToMerge); + + public: + int m_numVerts; + int m_numElems; + + CellTopologyData m_topology; + CellMapData m_mapData; + }; +} +} + +#endif \ No newline at end of file diff --git a/src/axom/mir/CellGenerator.cpp b/src/axom/mir/CellGenerator.cpp new file mode 100644 index 0000000000..3e5537325f --- /dev/null +++ b/src/axom/mir/CellGenerator.cpp @@ -0,0 +1,220 @@ +// 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 "CellGenerator.hpp" + + +namespace axom +{ +namespace mir +{ + +//-------------------------------------------------------------------------------- + +CellGenerator::CellGenerator() +{ + +} + +//-------------------------------------------------------------------------------- + +CellGenerator::~CellGenerator() +{ + +} + +//-------------------------------------------------------------------------------- + +void CellGenerator::generateTopologyData(const std::map >& newElements, + const std::map >& newVertices, + CellData& out_cellData) +{ + // Store the evInds and evBegins data in the output vectors + int currentEVBeginIndex = 0; + for (auto itr = newElements.begin(); itr != newElements.end(); itr++) + { + // Push the start index of the next element + out_cellData.m_topology.m_evBegins.push_back(currentEVBeginIndex); + + // Push the next element's vertices + for (unsigned int vIndex = 0; vIndex < itr->second.size(); ++vIndex) + { + out_cellData.m_topology.m_evInds.push_back(itr->second[vIndex]); + ++currentEVBeginIndex; + } + } + + // Push the index that occurs after the last vertex + out_cellData.m_topology.m_evBegins.push_back(currentEVBeginIndex); + + // Store the veInds and veBegins data in the output vectors + int currentVEBeginIndex = 0; + for (auto itr = newVertices.begin(); itr != newVertices.end(); itr++) + { + // Push the start index of the vertex's elements + out_cellData.m_topology.m_veBegins.push_back(currentVEBeginIndex); + + // Push the next vertex's elements + for (unsigned int eIndex = 0; eIndex < itr->second.size(); eIndex++) + { + out_cellData.m_topology.m_veInds.push_back(itr->second[eIndex]); + ++currentVEBeginIndex; + } + } + + // Push the index that occurs after the last element + out_cellData.m_topology.m_veBegins.push_back(currentVEBeginIndex); +} + +//-------------------------------------------------------------------------------- + +void CellGenerator::generateVertexPositions(const mir::Shape shapeType, + const std::map>& newVertices, + const std::vector& vertexPositions, + axom::float64* tValues, + CellData& out_cellData) +{ + for (auto itr = newVertices.begin(); itr != newVertices.end(); itr++) + { + int vID = itr->first; + + if ( vID < mir::utilities::numVerts(shapeType) ) + { + // This vertex is one of the shape's original vertices + out_cellData.m_mapData.m_vertexPositions.push_back( vertexPositions[vID] ); + } + else + { + // This vertex is between two of the shape's original vertices + int vIDFrom = mir::utilities::getEdgeEndpoint(shapeType, vID, true); + int vIDTo = mir::utilities::getEdgeEndpoint(shapeType, vID, false); + + out_cellData.m_mapData.m_vertexPositions.push_back( + mir::Point2::lerp( vertexPositions[vIDFrom], vertexPositions[vIDTo], tValues[vID] ) ); + } + } +} + +//-------------------------------------------------------------------------------- + +void CellGenerator::generateVertexVolumeFractions(const mir::Shape shapeType, + const std::map>& newVertices, + const std::vector >& vertexVF, + axom::float64* tValues, + CellData& out_cellData) +{ + out_cellData.m_mapData.m_vertexVolumeFractions.resize(vertexVF.size()); // vertexVF size is the number of materials + + for (auto itr = newVertices.begin(); itr != newVertices.end(); itr++) + { + int vID = itr->first; + + for (unsigned long matID = 0; matID < vertexVF.size(); ++matID) + { + if ( vID < mir::utilities::numVerts(shapeType) ) + { + // This vertex is one of the shape's original vertices + out_cellData.m_mapData.m_vertexVolumeFractions[matID].push_back( vertexVF[matID][vID] ); + } + else + { + // This vertex is between two of the shape's original vertices + int vIDFrom = mir::utilities::getEdgeEndpoint(shapeType, vID, true); + int vIDTo = mir::utilities::getEdgeEndpoint(shapeType, vID, false); + + out_cellData.m_mapData.m_vertexVolumeFractions[matID].push_back( + axom::utilities::lerp( vertexVF[matID][vIDFrom], vertexVF[matID][vIDTo], tValues[vID] ) ); + } + } + } +} + +//-------------------------------------------------------------------------------- + +int CellGenerator::determineCleanCellMaterial(const Shape elementShape, + const std::vector& vertexIDs, + const int matOne, + const int matTwo, + const std::vector >& vertexVF) +{ + int dominantMaterial = matOne; + + axom::float64 matOneVF = -1.0; + axom::float64 matTwoVF = -1.0; + + for (unsigned long it = 0; it < vertexIDs.size(); ++it) + { + int vID = vertexIDs[it]; + + if ( vID >= 0 && vID < mir::utilities::numVerts(elementShape) ) + { + if (matOne != NULL_MAT) + { + matOneVF = vertexVF[matOne][vID]; + } + if (matTwo != NULL_MAT) + { + matTwoVF = vertexVF[matTwo][vID]; + } + + dominantMaterial = (matOneVF > matTwoVF) ? matOne : matTwo; + } + } + + return dominantMaterial; +} + +//-------------------------------------------------------------------------------- + +mir::Shape CellGenerator::determineElementShapeType(const Shape parentShapeType, + const int numVerts) +{ + mir::Shape newShapeType; + if (parentShapeType == mir::Shape::Triangle || parentShapeType == mir::Shape::Quad) + { + // Handle the two-dimensional case + switch (numVerts) + { + case 3: + newShapeType = mir::Shape::Triangle; + break; + case 4: + newShapeType = mir::Shape::Quad; + break; + default: + newShapeType = mir::Shape::Triangle; + printf("Invalid number of vertices in determineElementShapeType().\n"); + break; + } + } + else + { + // Handle the three-dimensional case + switch (numVerts) + { + case 4: + newShapeType = mir::Shape::Tetrahedron; + break; + case 5: + newShapeType = mir::Shape::Pyramid; + break; + case 6: + newShapeType = mir::Shape::Triangular_Prism; + break; + case 8: + newShapeType = mir::Shape::Hexahedron; + break; + default: + newShapeType = mir::Shape::Tetrahedron; + printf("Invalid number of vertices in determineElementShapeType().\n"); + break; + } + } + return newShapeType; +} + + +} +} diff --git a/src/axom/mir/CellGenerator.hpp b/src/axom/mir/CellGenerator.hpp new file mode 100644 index 0000000000..ccf189537b --- /dev/null +++ b/src/axom/mir/CellGenerator.hpp @@ -0,0 +1,150 @@ +// 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) + +/** + * \file CellGenerator.hpp + * + * \brief Contains the specification for the CellGenerator class. + * + */ + +#ifndef __CELL_GENERATOR_H__ +#define __CELL_GENERATOR_H__ + +#include "axom/core.hpp" // for axom macros +#include "axom/slam.hpp" // unified header for slam classes and functions + +#include "MIRMesh.hpp" +#include "MIRUtilities.hpp" +#include "MIRMeshTypes.hpp" +#include "CellData.hpp" +#include "ZooClippingTables.hpp" +#include "MIRUtilities.hpp" +#include "CellClipper.hpp" + +//-------------------------------------------------------------------------------- + +namespace axom +{ +namespace mir +{ + +//-------------------------------------------------------------------------------- + + /** + * \class CellGenerator + * + * \brief A class that generates that uses the clipping information to generate + * the data needed in order to produce a clean, reconstructed mesh. + */ + class CellGenerator + { + public: + + /** + * \brief Default constructor. + */ + CellGenerator(); + + /** + * \brief Default destructor. + */ + ~CellGenerator(); + + /** + * \brief Generates the topology of the new elements resulting from a split. + * + * \param newElements An ordered map of the generated elements' IDs to a list of the vertex IDs in the local frame of the output element. + * \param newVertices An order map of the vertex IDs in the local frame of the output element to the generated elements' IDs it is associated with. + * \param out_cellData Container to store the topology data of the generated elements. + */ + void generateTopologyData(const std::map >& newElements, + const std::map >& newVertices, + CellData& out_cellData); + + + /** + * \brief Generates the vertex positions values for each of the new vertices of the generated element. + * + * \param shapeType The shape type of the element. + * \param newVertices An ordered map of vertices that compose the newly generated elements. + * \param vertexPositions A vector of positions for each of the original element's vertices. + * \param tValues An array of t values where each of the midpoint vertices are for the newly generated elements. + * \param out_cellData Container to store the vertex position data of the generated elements. + * + * \note New vertex positions are interpolated from the original vertex positions. + */ + void generateVertexPositions(const mir::Shape shapeType, + const std::map>& newVertices, + const std::vector& vertexPositions, + axom::float64* tValues, + CellData& out_cellData); + + + /** + * \brief Generates the vertex volume fractions for each of the new vertices of the generated element. + * + * \param shapeType The shape type of the element. + * \param newVertices An ordered map of vertices that compose the newly generated elements. + * \param vertexVF A vector of positions for each of the original element's vertices. + * \param tValues An array of t values where each of the midpoint vertices are for the newly generated elements. + * \param out_cellData Container to store the vertex position data of the generated elements. + * + * \note New vertex positions are interpolated from the original vertex volume fractions. + */ + void generateVertexVolumeFractions(const mir::Shape shapeType, + const std::map>& newVertices, + const std::vector >& vertexVF, + axom::float64* tValues, + CellData& out_cellData); + + /** + * \brief Determines the more dominant material of the two given for the given element. + * + * \param shapeType An enumerator denoting the element's shape. + * \param vertexIDs A list of vertex IDs into the vertexVF param. + * \param matOne The ID of the first material. + * \param matTwo The ID of the second material. + * \param vertexVF The list of volume fractions associated with the given vertices in the vertexIDs param. + * + * \return The ID of the dominant material of the element. + * + * \note The dominant element for the 2D/3D cases will be the same as the material present at one of the + * original vertices that existed prior to the split. So, if you can find this vertex and its dominant + * material, then you know the dominant material of this new element. + * + * \note It is assumed that the given cell is one that results from splitting its parent cell. + */ + int determineCleanCellMaterial(const Shape shapeType, + const std::vector& vertexIDs, + const int matOne, + const int matTwo, + const std::vector >& vertexVF); + + /** + * \brief Determine the shape type of an element. + * + * \param parentShapeType The shape of the element from which the new element is generated. + * \param numVerts The number of vertices of the new element. + * + * \note It is assumed that the given cell is one that results from splitting its parent cell. + */ + mir::Shape determineElementShapeType(const Shape parentShapeType, + const int numVerts); + + /** + * \brief Ensures that the element is dominated by a material that is actually present in the original parent cell. + * + * \param + */ + // void fixDominantMaterial(mir::MIRMesh& originalMesh, const std::map >& newElements, const int matOne, const int matTwo, CellData& out_cellData); + }; + +//-------------------------------------------------------------------------------- + +} +} + +#endif \ No newline at end of file diff --git a/src/axom/mir/InterfaceReconstructor.cpp b/src/axom/mir/InterfaceReconstructor.cpp new file mode 100644 index 0000000000..b18c17cfdb --- /dev/null +++ b/src/axom/mir/InterfaceReconstructor.cpp @@ -0,0 +1,274 @@ +// 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 "InterfaceReconstructor.hpp" + +namespace axom +{ +namespace mir +{ + +//-------------------------------------------------------------------------------- + +InterfaceReconstructor::InterfaceReconstructor() +{ + +} + +//-------------------------------------------------------------------------------- + +InterfaceReconstructor::~InterfaceReconstructor() +{ + +} + +//-------------------------------------------------------------------------------- + +void InterfaceReconstructor::computeReconstructedInterface(mir::MIRMesh& inputMesh, + mir::MIRMesh& outputMesh) +{ + // Store a reference to the original mesh + m_originalMesh = inputMesh; + + // Initialize the final mesh to be the same as the input mesh + mir::MIRMesh finalMesh(m_originalMesh); + + // For each material in the mesh, split the mesh on the current material and generate a new mesh (input into next iteration) + for (int matID = 0; matID < m_originalMesh.m_numMaterials; ++matID) + { + // Copy the mesh to be split + mir::MIRMesh intermediateMesh(finalMesh); + + // Create an array to store the output of each element being split. + CellData temp_cellData[intermediateMesh.m_elems.size()]; + + // Process/split each element + for (int eID = 0; eID < intermediateMesh.m_elems.size(); ++eID) + { + // Update the materials upon which the split will occur + int currentDominantMat = intermediateMesh.m_elementDominantMaterials[eID]; + int matOne = matID; + + generateCleanCells(eID, currentDominantMat, matOne, intermediateMesh, temp_cellData[eID]); + } + + // Merge each of the cells into the first CellData struct + for (int eID = 1; eID < intermediateMesh.m_elems.size(); ++eID) + { + temp_cellData[0].mergeCell(temp_cellData[eID]); + } + + mir::VertSet combined_verts(temp_cellData[0].m_numVerts); + mir::ElemSet combined_elems(temp_cellData[0].m_numElems); + + // Create the final, processed mesh + mir::MIRMesh processedMesh; + processedMesh.initializeMesh(combined_verts, combined_elems, intermediateMesh.m_numMaterials, temp_cellData[0].m_topology, temp_cellData[0].m_mapData); + processedMesh.constructMeshVolumeFractionsVertex(temp_cellData[0].m_mapData.m_vertexVolumeFractions); + + // Store the current mesh to be passed into the next iteration + finalMesh = processedMesh; + finalMesh.constructMeshRelations(); + } + + // TODO: Run post-processing on mesh to meld the mesh vertices that are within a certain epsilon + outputMesh = finalMesh; + // return finalMesh; +} + +//-------------------------------------------------------------------------------- + +void InterfaceReconstructor::computeReconstructedInterfaceIterative(mir::MIRMesh& inputMesh, + const int numIterations, + const axom::float64 percent, + mir::MIRMesh& outputMesh) +{ + int numElems = inputMesh.m_elems.size(); + int numMaterials = inputMesh.m_numMaterials; + + // Make a copy of the original input mesh + mir::MIRMesh meshToImprove(inputMesh); + + // Calculate the reconstruction on the unmodified, input mesh + computeReconstructedInterface(meshToImprove, outputMesh); + + for (int it = 0; it < numIterations; ++it) + { + printf("iteration %d\n", it); + // Calculate the output element volume fractions of the resulting output mesh + std::vector > resultingElementVF = outputMesh.computeOriginalElementVolumeFractions(); + + // Initialize the vector to store the improved element volume fractions + std::vector > improvedElementVF; + improvedElementVF.resize(numMaterials); + + // Modify the copy of the original mesh's element volume fractions by percent difference (modify meshToImprove's elementVF) + for (int matID = 0; matID < numMaterials; ++matID) + { + for (int eID = 0; eID < numElems; ++eID) + { + axom::float64 difference = inputMesh.m_materialVolumeFractionsElement[matID][eID] - resultingElementVF[matID][eID]; + improvedElementVF[matID].push_back( inputMesh.m_materialVolumeFractionsElement[matID][eID] + ( difference * percent ) ); + } + } + + // Reconstruct the element AND vertex VFs + meshToImprove.constructMeshVolumeFractionsMaps(improvedElementVF); + + // Calculate the reconstruction on the modified, input mesh + computeReconstructedInterface(meshToImprove, outputMesh); + } +} + +//-------------------------------------------------------------------------------- + +void InterfaceReconstructor::generateCleanCells(const int eID, + const int matOne, + const int matTwo, + mir::MIRMesh& tempMesh, + CellData& out_cellData) +{ + mir::Shape shapeType = (mir::Shape) tempMesh.m_shapeTypes[eID]; + auto elementVertices = tempMesh.m_bdry[eID]; + + // Set up data structures needed to compute how element should be clipped + std::map > newElements; // hashmap of the new elements' vertices | Note: maps are ordered sets + std::map > newVertices; // hashmap of the new vertices' elements | Note: maps are ordered sets + std::vector verticesPresent(mir::utilities::maxPossibleNumVerts(shapeType), 0); // Vector of flags denoting whether the vertex is present in the current case or not + axom::float64* tValues = new axom::float64[ mir::utilities::maxPossibleNumVerts(shapeType) ]{0}; // Array of t values that denote the percent value of where the edge should be clipped + + + // Set up the volume fractions for the current element for the two materials currently being considered + std::vector > vertexVF(2); + for (int vID = 0; vID < elementVertices.size(); ++vID) + { + int originalVID = elementVertices[vID]; + if (matOne == NULL_MAT) + { + vertexVF[0].push_back(-1.0); + } + else + { + vertexVF[0].push_back( tempMesh.m_materialVolumeFractionsVertex[matOne][originalVID] ); + } + + if (matTwo == NULL_MAT) + { + vertexVF[1].push_back(-1.0); + } + else + { + vertexVF[1].push_back( tempMesh.m_materialVolumeFractionsVertex[matTwo][originalVID] ); + } + } + + // Clip the element + CellClipper clipper; + clipper.computeClippingPoints(shapeType, vertexVF, newElements, newVertices, tValues); + + // Determine which vertices are present + for (auto itr = newVertices.begin(); itr != newVertices.end(); itr++) + { + int vID = itr->first; + verticesPresent[vID] = 1; + } + + // Calculate the total number of elements and vertices that were generated from splitting the current element + out_cellData.m_numElems = (int) newElements.size(); + out_cellData.m_numVerts = (int) newVertices.size(); + + // Generate the topology and connectivity of the newly split element + CellGenerator cellGenerator; + cellGenerator.generateTopologyData(newElements, newVertices, out_cellData); + + // Generate the vertex position values of the newly split element + std::vector originalElementVertexPositions; + for (int vID = 0; vID < elementVertices.size(); ++vID) + { + int originalVID = elementVertices[vID]; + originalElementVertexPositions.push_back( tempMesh.m_vertexPositions[originalVID] ); + } + cellGenerator.generateVertexPositions( shapeType, newVertices, originalElementVertexPositions, tValues, out_cellData); + + // Generate the vertex volume fractions of the newly split elements + std::vector > originalElementVertexVF; + for (int matID = 0; matID < tempMesh.m_numMaterials; ++matID) + { + std::vector materialVertexVF; + for (int vID = 0; vID < elementVertices.size(); ++vID) + { + int originalVID = elementVertices[vID]; + + materialVertexVF.push_back( tempMesh.m_materialVolumeFractionsVertex[matID][originalVID] ); + } + originalElementVertexVF.push_back( materialVertexVF ); + } + + cellGenerator.generateVertexVolumeFractions( shapeType, newVertices, originalElementVertexVF, tValues, out_cellData); + + // Determine and store the dominant material of this element + for (auto itr = newElements.begin(); itr != newElements.end(); itr++) + { + int dominantMaterial = cellGenerator.determineCleanCellMaterial(shapeType, itr->second, matOne, matTwo, out_cellData.m_mapData.m_vertexVolumeFractions); + out_cellData.m_mapData.m_elementDominantMaterials.push_back(dominantMaterial); + } + + // Determine and store the parent of this element + out_cellData.m_mapData.m_elementParents.resize(newElements.size()); + for (auto itr = newElements.begin(); itr != newElements.end(); itr++) + { + out_cellData.m_mapData.m_elementParents[itr->first] = tempMesh.m_elementParentIDs[eID]; + } + + // Determine the generated elements' shape types + out_cellData.m_mapData.m_shapeTypes.resize(newElements.size()); + for (auto itr = newElements.begin(); itr != newElements.end(); itr++) + { + out_cellData.m_mapData.m_shapeTypes[itr->first] = cellGenerator.determineElementShapeType(shapeType, itr->second.size()); + } + + // TODO: Make a function that does this + // Check that each element is dominated by a material that is actually present in the original parent cell + for (auto itr = newElements.begin(); itr != newElements.end(); itr++) + { + int parentElementID = out_cellData.m_mapData.m_elementParents[itr->first]; + int currentDominantMaterial = out_cellData.m_mapData.m_elementDominantMaterials[itr->first]; + + if (m_originalMesh.m_materialVolumeFractionsElement[currentDominantMaterial][parentElementID] == 0) + { + // This material is not present in the original element from which the current element comes from + if (currentDominantMaterial == matOne) + out_cellData.m_mapData.m_elementDominantMaterials[itr->first] = matTwo; + else + out_cellData.m_mapData.m_elementDominantMaterials[itr->first] = matOne; + } + } + + // Modify the evIndex values to account for the fact that perhaps not all 8 possible vertices are present + std::vector evIndexSubtract; + evIndexSubtract.resize(mir::utilities::maxPossibleNumVerts(shapeType), 0); + for (unsigned long i = 0; i < evIndexSubtract.size(); ++i) + { + if (verticesPresent[i] == 1) + evIndexSubtract[i] = 0; + else + evIndexSubtract[i] = 1; + } + for (unsigned long i = 1; i < evIndexSubtract.size(); ++i) + evIndexSubtract[i] += evIndexSubtract[i - 1]; + + for (unsigned long i = 0; i < out_cellData.m_topology.m_evInds.size(); ++i) + { + out_cellData.m_topology.m_evInds[i] -= evIndexSubtract[ out_cellData.m_topology.m_evInds[i] ]; + } + + // Memory management + delete[] tValues; +} + +//-------------------------------------------------------------------------------- + +} +} diff --git a/src/axom/mir/InterfaceReconstructor.hpp b/src/axom/mir/InterfaceReconstructor.hpp new file mode 100644 index 0000000000..ed674af4a4 --- /dev/null +++ b/src/axom/mir/InterfaceReconstructor.hpp @@ -0,0 +1,107 @@ +// 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) + +/** + * \file InterfaceReconstructor.hpp + * + * \brief Contains the specification for the InterfaceReconstructor class. + * + */ + +#ifndef __INTERFACE_RECONSTRUCTOR_H__ +#define __INTERFACE_RECONSTRUCTOR_H__ + +#include "axom/core.hpp" // for axom macros +#include "axom/slam.hpp" + +#include "MIRMesh.hpp" +#include "CellData.hpp" +#include "ZooClippingTables.hpp" +#include "MIRUtilities.hpp" +#include "CellClipper.hpp" +#include "CellGenerator.hpp" + +#include + +namespace numerics = axom::numerics; +namespace slam = axom::slam; + +namespace axom +{ +namespace mir +{ + + /** + * \class InterfaceReconstructor + * + * \brief A class that contains the functionality for taking an input mesh + * with mixed cells and outputs a mesh with clean cells. + * + * \detail This class requires that the user create a mesh of type MIRMesh + * and pass that into one of the reconstruction methods. There are + * currently two reconstructions methods implemented: one is the + * a zoo-based algorithm described in Meredith 2004, and the other + * is the iterative version described in Meredith and Childs 2010. + */ + class InterfaceReconstructor + { + public: + + /** + * \brief Default constructor. + */ + InterfaceReconstructor(); + + + /** + * \brief Default destructor. + */ + ~InterfaceReconstructor(); + + + /** + * \brief Performs material interface reconstruction using the zoo-based algorithm. + * + * \param inputMesh The mesh composed of mixed cells. + * \param outputMesh The mesh composed of clean cells. + */ + void computeReconstructedInterface(mir::MIRMesh& inputMesh, + mir::MIRMesh& outputMesh); + + + /** + * \brief Performs material interface reconstruction using an iterative version of the zoo-based algorithm. + * + * \param inputMesh The mesh made up of mixed cells. + * \param numIterations The number of iterations for which to run the algorithm. + * \param percent The percent of the difference to use when modifying the original element volume fractions. + * \param outputMesh The mesh composed of clean cells. + */ + void computeReconstructedInterfaceIterative(mir::MIRMesh& inputMesh, + const int numIterations, + const axom::float64 percent, + mir::MIRMesh& outputMesh); + + /** + * \brief Generates a set of clean cells by splitting the element with the two given materials. + * + * \param eID The ID of the element to be split. + * \param matOne The first material to split with. + * \param matTwo The second material to split with. + * \param tempMesh A reference to the mesh the element comes from. + * \param out_cellData Container to store the data of the generated elements. + */ + void generateCleanCells(const int eID, + const int matOne, + const int matTwo, + mir::MIRMesh& tempMesh, + CellData& out_cellData); + + private: + mir::MIRMesh m_originalMesh; + }; +} +} +#endif \ No newline at end of file diff --git a/src/axom/mir/MIRMesh.cpp b/src/axom/mir/MIRMesh.cpp new file mode 100644 index 0000000000..9181c10522 --- /dev/null +++ b/src/axom/mir/MIRMesh.cpp @@ -0,0 +1,674 @@ +// 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 "MIRMesh.hpp" + +#include "axom/core.hpp" +#include "axom/primal.hpp" + +namespace axom +{ +namespace mir +{ + + +MIRMesh::MIRMesh() + : m_verts(0) + , m_elems(0) + , m_numMaterials(0) +{} + +//-------------------------------------------------------------------------------- + +MIRMesh::MIRMesh(const MIRMesh& other) + : m_verts(other.m_verts) + , m_elems(other.m_elems) + , m_materialVolumeFractionsElement(other.m_materialVolumeFractionsElement) + , m_materialVolumeFractionsVertex(other.m_materialVolumeFractionsVertex) + , m_numMaterials(other.m_numMaterials) + , m_meshTopology(other.m_meshTopology) +{ + constructMeshRelations(); + constructVertexPositionMap(other.m_vertexPositions.data()); + constructElementParentMap(other.m_elementParentIDs.data()); + constructElementDominantMaterialMap(other.m_elementDominantMaterials.data()); + + m_shapeTypes= IntMap( &m_elems ); + for(auto el: m_elems) { m_shapeTypes[el] = other.m_shapeTypes[el]; } + +} + +MIRMesh& MIRMesh::operator=(const MIRMesh& other) +{ + if(this != &other) + { + m_verts = other.m_verts; + m_elems = other.m_elems; + + m_meshTopology = other.m_meshTopology; + + constructMeshRelations(); + constructVertexPositionMap(other.m_vertexPositions.data()); + constructElementParentMap(other.m_elementParentIDs.data()); + constructElementDominantMaterialMap(other.m_elementDominantMaterials.data()); + + + m_shapeTypes = IntMap( &m_elems ); + for(auto el: m_elems) { m_shapeTypes[el] = other.m_shapeTypes[el]; } + + m_materialVolumeFractionsElement = other.m_materialVolumeFractionsElement; + m_materialVolumeFractionsVertex = other.m_materialVolumeFractionsVertex; + + m_numMaterials = other.m_numMaterials; + } + return *this; +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::initializeMesh(const VertSet _verts, + const ElemSet _elems, + const int _numMaterials, + const CellTopologyData& _topology, + const CellMapData& _mapData, + const std::vector >& _elementVF) +{ + // Initialize the vertex and element sets + m_verts = _verts; + m_elems = _elems; + + m_numMaterials = _numMaterials; + + // Intialize the mesh topology + m_meshTopology.m_evInds = _topology.m_evInds; + m_meshTopology.m_evBegins = _topology.m_evBegins; + m_meshTopology.m_veInds = _topology.m_veInds; + m_meshTopology.m_veBegins = _topology.m_veBegins; + + // Initialize the mesh relations + constructMeshRelations(); + + // Initialize the mesh's data maps + constructVertexPositionMap(_mapData.m_vertexPositions); + constructElementParentMap(_mapData.m_elementParents); + constructElementDominantMaterialMap(_mapData.m_elementDominantMaterials); + constructElementShapeTypesMap(_mapData.m_shapeTypes); + + // Initialize the element and vertex volume fraction maps + if (_elementVF.size() > 0) + constructMeshVolumeFractionsMaps(_elementVF); + + // Check validity of the sets + SLIC_ASSERT_MSG( m_verts.isValid(), "Vertex set is not valid."); + SLIC_ASSERT_MSG( m_elems.isValid(), "Element set is not valid."); +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::constructMeshRelations() +{ + // construct boundary relation from elements to vertices using variable cardinality + { + using RelationBuilder = ElemToVertRelation::RelationBuilder; + m_bdry = RelationBuilder() + .fromSet( &m_elems ) + .toSet( &m_verts ) + .begins( RelationBuilder::BeginsSetBuilder() + .size( m_elems.size() ) + .data( m_meshTopology.m_evBegins.data() ) ) + .indices ( RelationBuilder::IndicesSetBuilder() + .size( m_meshTopology.m_evInds.size() ) + .data( m_meshTopology.m_evInds.data() ) ); + } + + + { + // _quadmesh_example_construct_cobdry_relation_start + // construct coboundary relation from vertices to elements + using RelationBuilder = VertToElemRelation::RelationBuilder; + m_cobdry = RelationBuilder() + .fromSet( &m_verts ) + .toSet( &m_elems ) + .begins( RelationBuilder::BeginsSetBuilder() + .size( m_verts.size() ) + .data( m_meshTopology.m_veBegins.data() ) ) + .indices( RelationBuilder::IndicesSetBuilder() + .size( m_meshTopology.m_veInds.size() ) + .data( m_meshTopology.m_veInds.data() ) ); + // _quadmesh_example_construct_cobdry_relation_end + } + + SLIC_ASSERT_MSG( m_bdry.isValid(), "Boundary relation is not valid."); + SLIC_ASSERT_MSG( m_cobdry.isValid(), "Coboundary relation is not valid."); + + SLIC_DEBUG("Elem-Vert relation has size " << m_bdry.totalSize()); + SLIC_DEBUG("Vert-Elem relation has size " << m_cobdry.totalSize()); +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::constructMeshVolumeFractionsMaps(const std::vector >& elementVF) +{ + // Clear the old maps + m_materialVolumeFractionsElement.clear(); + m_materialVolumeFractionsVertex.clear(); + + // Initialize the maps for all of the materials with the input volume fraction data for each material + for (int matID = 0; matID < m_numMaterials; ++matID) + { + // Initialize the map for the current material + m_materialVolumeFractionsElement.push_back(ScalarMap( &m_elems )); + + // Copy the data for the current material + for (int eID = 0; eID < m_elems.size(); ++eID) + { + m_materialVolumeFractionsElement[matID][eID] = elementVF[matID][eID]; + } + + SLIC_ASSERT_MSG( m_materialVolumeFractionsElement[matID].isValid(), "Element volume fraction map is not valid."); + } + + // Initialize the maps for all of the vertex volume fractions + for (int matID = 0; matID < m_numMaterials; ++matID) + { + // Initialize the new map for the volume fractions + m_materialVolumeFractionsVertex.push_back(ScalarMap( &m_verts ) ); + + // Calculate the average volume fraction value for the current vertex for the current material + for (int vID = 0; vID < m_verts.size(); ++vID) + { + // Compute the per vertex volume fractions for the green material + axom::float64 sum = 0; + auto vertexElements = m_cobdry[vID]; + + for (int i = 0; i < vertexElements.size(); ++i) + { + auto eID = vertexElements[i]; + sum += m_materialVolumeFractionsElement[matID][eID]; + } + + m_materialVolumeFractionsVertex[matID][vID] = sum / vertexElements.size(); + } + + SLIC_ASSERT_MSG( m_materialVolumeFractionsVertex[matID].isValid(), "Vertex volume fraction map is not valid."); + } +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::constructMeshVolumeFractionsVertex(const std::vector >& vertexVF) +{ + // Initialize the maps for all of the materials with the input volume fraction data for each vertex + for (int matID = 0; matID < m_numMaterials; ++matID) + { + // Initialize the map for the current material + m_materialVolumeFractionsVertex.push_back(ScalarMap( &m_verts )); + + // Copy the data for the current material + for (int vID = 0; vID < m_verts.size(); ++vID) + { + m_materialVolumeFractionsVertex[matID][vID] = vertexVF[matID][vID]; + } + + SLIC_ASSERT_MSG( m_materialVolumeFractionsVertex[matID].isValid(), "Vertex volume fraction map is not valid."); + } +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::constructVertexPositionMap(const std::vector& data) +{ + // construct the position map on the vertices + m_vertexPositions = PointMap( &m_verts ); + + for (int vID = 0; vID < m_verts.size(); ++vID) + m_vertexPositions[vID] = data[vID]; + + SLIC_ASSERT_MSG( m_vertexPositions.isValid(), "Position map is not valid."); +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::constructElementParentMap(const std::vector& elementParents) +{ + // Initialize the map for the elements' parent IDs + m_elementParentIDs = IntMap( &m_elems ); + + // Copy the data for the elements + for (int eID = 0; eID < m_elems.size(); ++eID) + m_elementParentIDs[eID] = elementParents[eID]; + + SLIC_ASSERT_MSG( m_elementParentIDs.isValid(), "Element parent map is not valid."); +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::constructElementDominantMaterialMap(const std::vector& dominantMaterials) +{ + // Initialize the map for the elements' dominant colors + m_elementDominantMaterials = IntMap( &m_elems ); + + // Copy the dat for the elements + for (int eID = 0; eID < m_elems.size(); ++eID) + m_elementDominantMaterials[eID] = dominantMaterials[eID]; + + SLIC_ASSERT_MSG( m_elementDominantMaterials.isValid(), "Element dominant materials map is not valid."); +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::constructElementShapeTypesMap(const std::vector& shapeTypes) +{ + // Initialize the map for the elements' dominant colors + m_shapeTypes = IntMap( &m_elems ); + + // Copy the data for the elements + for (int eID = 0; eID < m_elems.size(); ++eID) + m_shapeTypes[eID] = shapeTypes[eID]; + + SLIC_ASSERT_MSG( m_shapeTypes.isValid(), "Element dominant materials map is not valid."); +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::print() +{ + printf("\n------------------------Printing Mesh Information:------------------------\n"); + printf("number of vertices: %d\n", m_verts.size()); + printf("number of elements: %d\n", m_elems.size()); + printf("number of materials: %d\n", m_numMaterials); + + printf("evInds: { "); + for (unsigned long i = 0; i < m_meshTopology.m_evInds.size(); i++) + { + printf("%d ", m_meshTopology.m_evInds[i]); + } + printf("}\n"); + + printf("evBegins: { "); + for (unsigned long i = 0; i < m_meshTopology.m_evBegins.size(); i++) + { + printf("%d ", m_meshTopology.m_evBegins[i]); + } + printf("}\n"); + + printf("veInds: { "); + for (unsigned long i = 0; i < m_meshTopology.m_veInds.size(); i++) + { + printf("%d ", m_meshTopology.m_veInds[i]); + } + printf("}\n"); + + printf("veBegins: { "); + for (unsigned long i = 0; i < m_meshTopology.m_veBegins.size(); i++) + { + printf("%d ", m_meshTopology.m_veBegins[i]); + } + printf("}\n"); + + printf("vertexPositions: { "); + for (int i = 0; i < m_verts.size(); ++i) + { + printf("{%.2f, %.2f} ", m_vertexPositions[i][0], m_vertexPositions[i][1]); + } + printf("}\n"); + + printf("elementParentIDs: { "); + for (int i = 0; i < m_elems.size(); ++i) + { + printf("%d ", m_elementParentIDs[i]); + } + printf("}\n"); + + printf("elementDominantMaterials: { "); + for (int i = 0; i < m_elems.size(); ++i) + { + printf("%d ", m_elementDominantMaterials[i]); + } + printf("}\n"); + + printf("shapeTypes: { "); + for (int i = 0; i < m_elems.size(); ++i) + { + printf("%d ", m_shapeTypes[i]); + } + printf("}\n"); + + printf("elementVolumeFractions: { \n"); + for (unsigned long i = 0; i < m_materialVolumeFractionsElement.size(); ++i) + { + printf(" { "); + for (int j = 0; j < m_elems.size(); ++j) + { + printf("%.3f, ", m_materialVolumeFractionsElement[i][j]); + } + printf("}\n"); + } + printf("}\n"); + + printf("vertexVolumeFractions: { \n"); + for (unsigned long i = 0; i < m_materialVolumeFractionsVertex.size(); ++i) + { + printf(" { "); + for (int j = 0; j < m_verts.size(); ++j) + { + printf("%.3f, ", m_materialVolumeFractionsVertex[i][j]); + } + printf("}\n"); + } + printf("}\n"); + printf("--------------------------------------------------------------------------\n"); +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::readMeshFromFile(std::string filename) +{ + printf("Mesh reading functionality not implemented yet. Can't read file: %s", filename.c_str()); + + // Read in header + + // Read in POINTS + + // Read in CELLS + + // Read in CELL_TYPES + + // Read in CELL_DATA (element volume fractions) +} + +//-------------------------------------------------------------------------------- + +void MIRMesh::writeMeshToFile(std::string filename) +{ + std::ofstream meshfile; + meshfile.open(filename); + std::ostream_iterator out_it(meshfile, " "); + + // write header + meshfile << "# vtk DataFile Version 3.0\n" + << "vtk output\n" + << "ASCII\n" + << "DATASET UNSTRUCTURED_GRID\n" + << "POINTS " << m_verts.size() << " double\n"; + + // write positions + for (int vID = 0; vID < m_verts.size(); ++vID) + { + meshfile << m_vertexPositions[vID][0] << " " << m_vertexPositions[vID][1] << " 0\n"; // must always set all 3 coords; set Z=0 for 2D + } + + meshfile << "\nCELLS " << m_elems.size() << " " << m_meshTopology.m_evInds.size() + m_elems.size(); + for (int i = 0; i < m_elems.size(); ++i) + { + int nVerts = m_meshTopology.m_evBegins[i+1] - m_meshTopology.m_evBegins[i]; + meshfile << "\n" << nVerts; + for (int j = 0; j < nVerts; ++j) + { + int startIndex = m_meshTopology.m_evBegins[i]; + meshfile << " " << m_meshTopology.m_evInds[startIndex + j]; + } + } + + meshfile << "\n\nCELL_TYPES " << m_elems.size() << "\n"; + for (int i = 0; i < m_elems.size(); ++i) + { + if (m_shapeTypes[i] == mir::Shape::Triangle) + meshfile << "5\n"; + else if (m_shapeTypes[i] == mir::Shape::Quad) + meshfile << "9\n"; + } + + // write element materials + meshfile << "\n\nCELL_DATA " << m_elems.size() + << "\nSCALARS materialIDs int 1" + << "\nLOOKUP_TABLE default \n"; + for(int i=0 ; i< m_elems.size() ; ++i) + { + meshfile << m_elementDominantMaterials[i] << " "; + } + + meshfile <<"\n"; +} + +//-------------------------------------------------------------------------------- + +std::vector > MIRMesh::computeOriginalElementVolumeFractions() +{ + std::map totalAreaOriginalElements; // the total area of the original elements + std::map newElementAreas; // the area of each of the generated child elements + std::map numChildren; // the number of child elements generated from one of the original elements + + // Compute the total area of each element of the original mesh and also the area of each new element + for (int eID = 0; eID < m_elems.size(); ++eID) + { + int numVertices = m_meshTopology.m_evBegins[eID + 1] - m_meshTopology.m_evBegins[eID]; + if (numVertices == 3) + { + Point2 trianglePoints[3]; + for (int i = 0; i < 3; ++i) + trianglePoints[i] = m_vertexPositions[m_meshTopology.m_evInds[ m_meshTopology.m_evBegins[eID] + i] ]; + newElementAreas[eID] = computeTriangleArea(trianglePoints[0], trianglePoints[1], trianglePoints[2]); + } + if (numVertices == 4) + { + Point2 trianglePoints[4]; + for (int i = 0; i < 4; ++i) + trianglePoints[i] = m_vertexPositions[m_meshTopology.m_evInds[ m_meshTopology.m_evBegins[eID] + i] ]; + newElementAreas[eID] = computeQuadArea(trianglePoints[0], trianglePoints[1], trianglePoints[2], trianglePoints[3]); + } + + totalAreaOriginalElements[ m_elementParentIDs[eID] ] += newElementAreas[eID]; + numChildren[ m_elementParentIDs[eID] ] += .4; + } + + // Intialize the element volume fraction vectors + std::vector > elementVolumeFractions; // indexed as: elementVolumeFractions[material][originalElementID] = volumeFraction + elementVolumeFractions.resize( m_numMaterials ); + for (unsigned long i = 0; i < elementVolumeFractions.size(); ++i) + elementVolumeFractions[i].resize( totalAreaOriginalElements.size(), 0.0 ); + + // Compute the volume fractions for each of the original mesh elements + for (auto itr = newElementAreas.begin(); itr != newElementAreas.end(); itr++) + { + int parentElementID = m_elementParentIDs[itr->first]; + int materialID = m_elementDominantMaterials[itr->first]; + elementVolumeFractions[materialID][parentElementID] += (itr->second / totalAreaOriginalElements[parentElementID]); + } + + return elementVolumeFractions; +} + +//-------------------------------------------------------------------------------- + +axom::float64 MIRMesh::computeTriangleArea(Point2 p0, + Point2 p1, + Point2 p2) +{ + return primal::Triangle(p0,p1,p2).area(); +} + +//-------------------------------------------------------------------------------- + +axom::float64 MIRMesh::computeQuadArea(Point2 p0, + Point2 p1, + Point2 p2, + Point2 p3) +{ + return computeTriangleArea(p0, p1, p2) + computeTriangleArea(p2, p3, p0); +} + +//-------------------------------------------------------------------------------- + +bool MIRMesh::isValid(bool verbose) const +{ + bool bValid = true; + + // Check the topology data + bValid = bValid && isTopologyValid(verbose); + + // Check the volume fraction data + bValid = bValid && areVolumeFractionsValid(verbose); + + // Check the map data + bValid = bValid && areMapsValid(verbose); + + return bValid; +} + +bool MIRMesh::isTopologyValid(bool verbose) const +{ + bool bValid = true; + + // check the vertex and element sets + bValid = bValid && m_verts.isValid(verbose); + bValid = bValid && m_elems.isValid(verbose); + + // check the relations + if(!m_verts.empty() && !m_elems.empty() ) + { + bValid = bValid && m_bdry.isValid(verbose); + bValid = bValid && m_cobdry.isValid(verbose); + } + + if(!bValid && verbose) + { + SLIC_WARNING("MIRMesh invalid: Invalid topology"); + } + + return bValid; +} + +bool MIRMesh::areVolumeFractionsValid(bool verbose) const +{ + bool bValid = true; + + int vfElemSize = m_materialVolumeFractionsElement.size(); + if( vfElemSize != m_numMaterials) + { + bValid = false; + if(verbose) + { + SLIC_WARNING("MIRMesh invalid: Invalid number of materials in volume fraction array."); + } + return bValid; + } + + // Check the sizes of the volume fraction arrays + for(int i=0; i < m_numMaterials; ++i) + { + const auto& mats = m_materialVolumeFractionsElement[i]; + if( mats.size() != m_elems.size()) + { + bValid = false; + if(verbose) + { + SLIC_WARNING("MIRMesh invalid: Material " << i <<" has wrong " + << " number of materials in volume fraction array." + << " Expected " << m_elems.size() + << ", got " << mats.size() << "."); + } + } + } + + if(!bValid) + { + return false; + } + + // Check that the volume fractions sum to 1.0 + for(auto el : m_elems.positions()) + { + double sum = 0; + for(int m=0; m < m_numMaterials; ++m) + { + sum += m_materialVolumeFractionsElement[m][el]; + } + + if(! axom::utilities::isNearlyEqual(sum, 1.)) + { + bValid = false; + if(verbose) + { + std::stringstream sstr; + + for(int m=0; m < m_numMaterials; ++m) + { + sstr << m_materialVolumeFractionsElement[m][el] << " "; + } + + SLIC_WARNING("MIRMesh invalid: Sum of materials in element" + << el << " was " << sum << ". Expected 1." + << " Values were: " << sstr.str() ); + } + } + } + + return bValid; +} + +bool MIRMesh::areMapsValid(bool verbose) const +{ + bool bValid = true; + + // Check the positions map + if( m_vertexPositions.size() != m_verts.size()) + { + bValid = false; + if(verbose) + { + SLIC_WARNING("MIRMesh invalid: Incorrect number of vertex positions." + << " Expected " << m_verts.size() + << ". Got " << m_vertexPositions.size()); + } + } + + bValid = bValid && m_vertexPositions.isValid(verbose); + + // Check the element maps + bValid = bValid && m_elementParentIDs.isValid(verbose); + bValid = bValid && m_elementDominantMaterials.isValid(verbose); + bValid = bValid && m_shapeTypes.isValid(verbose); + + if( m_elementParentIDs.size() != m_elems.size() ) + { + bValid = false; + if(verbose) + { + SLIC_WARNING("MIRMesh invalid: Incorrect number of elem parent IDs." + << " Expected " << m_elems.size() + << ". Got " << m_elementParentIDs.size()); + } + } + + if( m_elementDominantMaterials.size() != m_elems.size() ) + { + bValid = false; + if(verbose) + { + SLIC_WARNING("MIRMesh invalid: Incorrect number of elem dominant mats." + << " Expected " << m_elems.size() + << ". Got " << m_elementDominantMaterials.size()); + } + } + + if( m_shapeTypes.size() != m_elems.size() ) + { + bValid = false; + if(verbose) + { + SLIC_WARNING("MIRMesh invalid: Incorrect number of elem shape types." + << " Expected " << m_elems.size() + << ". Got " << m_shapeTypes.size()); + } + } + + return bValid; + +} + + +} +} diff --git a/src/axom/mir/MIRMesh.hpp b/src/axom/mir/MIRMesh.hpp new file mode 100644 index 0000000000..57b1e7446d --- /dev/null +++ b/src/axom/mir/MIRMesh.hpp @@ -0,0 +1,237 @@ +// 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) + +/** + * \file MIRMesh.hpp + * + * \brief Contains the specification for the MIRMesh class. + * + */ + +#ifndef __MIR_MESH_H__ +#define __MIR_MESH_H__ + + +#include "axom/core.hpp" // for axom macros +#include "axom/slam.hpp" // unified header for slam classes and functions + +#include "MIRMeshTypes.hpp" +#include "CellData.hpp" + +// C/C++ includes +#include // for definition of M_PI, exp() +#include +#include +#include + +namespace numerics = axom::numerics; +namespace slam = axom::slam; + +//-------------------------------------------------------------------------------- +namespace axom +{ +namespace mir +{ + + const int NULL_MAT = -1; + + //-------------------------------------------------------------------------------- + + /** + * \class MIRMesh + * + * \brief The MIRMesh class represents a finite element mesh containing element volume fractions. + * + * \detail This class is meant to be used in conjunction with the InterfaceReconstructor class + * to process the mesh. + * + */ + class MIRMesh + { + public: + /** + * \brief Default constructor. + */ + MIRMesh(); + + /** + * \brief Copy constructor. + */ + MIRMesh(const MIRMesh& other); + + /** + * \brief Default destructor. + */ + ~MIRMesh() = default; + + /** + * \brief Copy assignment. + */ + MIRMesh& operator=(const MIRMesh& other); + + /** + * \brief Initializes a mesh with the provided data and topology. + * + * \param _verts The set of vertices of the mesh. + * \param _elems The set of elements of the mesh. + * \param _numMaterials The number of materials present in the mesh. + * \param _topology The topology/connectivity of the mesh. + * \param _mapData The data used to initialized the maps associated with the vertex and element sets. + * \param _elementVF The volume fractions of each element. Note that this is an optional parameter. + */ + void initializeMesh(const VertSet _verts, + const ElemSet _elems, + const int _numMaterials, + const CellTopologyData& _topology, + const CellMapData& _mapData, + const std::vector >& _elementVF = {}); + + /** + * \brief Constructs the mesh boundary and coboundary relations. + * + * \note The boundary relation is from each element to its vertices. + * \note The coboundary relation is from each vertex to its elements. + */ + void constructMeshRelations(); + + /** + * \brief Constructs the element and vertex volume fraction maps. + * + * \param elementVF The volume fractions of each element. + */ + void constructMeshVolumeFractionsMaps(const std::vector >& elementVF); + + /** + * \brief Constructs the vertex volume fraction map. + * + * \param vertexVF The volume fractions of each vertex. + */ + void constructMeshVolumeFractionsVertex(const std::vector >& vertexVF); + + /** + * \brief Prints out the data contained within this mesh in a nice format. + */ + void print(); + + /** + * \brief Reads in a mesh specified by the given file. + * + * \param filename The location where the mesh file will be read from. + */ + void readMeshFromFile(const std::string filename); + + /** + * \brief Writes out the mesh to the given file. + * + * \param filename The location where the mesh file will be written. + * + * \note Currently reads in an ASCII, UNSTRUCTURED_GRID .vtk file. + */ + void writeMeshToFile(const std::string filename); + + + /** + * \brief Computes the volume fractions of the elements of the original mesh. + */ + std::vector > computeOriginalElementVolumeFractions(); + + /** + * \brief Checks that this instance of MIRMesh is valid + * + * \return True if this instance is valid, false otherwise + */ + bool isValid(bool verbose) const; + + private: + + /// Utility predicate to check if the mesh's topology is valid + bool isTopologyValid(bool verbose) const; + /// Utility predicate to check if the mesh's volume fractions are valid + bool areVolumeFractionsValid(bool verbose) const; + /// Utility predicate to check if the mesh's maps are valid + bool areMapsValid(bool verbose) const; + + /** + * \brief Constructs the positions map on the vertices. + * + * \param data The vector of position data for each vertex. + */ + void constructVertexPositionMap(const std::vector& data); + + /** + * \brief Constructs the map of elements to their original element parent. + * + * \param cellParents The vector of parent IDs for each element of the mesh. + */ + void constructElementParentMap(const std::vector& elementParents); + + /** + * \brief Constructs the map of elements to their dominant materials. + * + * \param dominantMaterials A vector of material ids that are the dominant material of each element. + */ + void constructElementDominantMaterialMap(const std::vector& dominantMaterials); + + /** + * \brief Constructs the map of elements to their shape types. + * + * \param shapeTypes A vector of shape enumerators that are the shape type of each element. + */ + void constructElementShapeTypesMap(const std::vector& shapeTypes); + + /** + * \brief Computes the area of the triangle defined by the given three vertex positions using Heron's formula. + * + * \param p0 The position of the first vertex. + * \param p1 The position of the second vertex. + * \param p2 The position of the third vertex. + */ + axom::float64 computeTriangleArea(Point2 p0, + Point2 p1, + Point2 p2); + + /** + * \brief Computes the area of the quad defined by the given four vertex positions. + * + * \param p0 The position of the first vertex. + * \param p1 The position of the second vertex. + * \param p2 The position of the third vertex. + * \param p3 The position of the fourth vertex. + * + * \note It is assumed that the points are given in consecutive, counter-clockwise order. + */ + axom::float64 computeQuadArea(Point2 p0, + Point2 p1, + Point2 p2, + Point2 p3); + + /**************************************************************** + * VARIABLES + ****************************************************************/ + public: + // Mesh Set Definitions + VertSet m_verts; // the set of vertices in the mesh + ElemSet m_elems; // the set of elements in the mesh + + // Mesh Relation Definitions + ElemToVertRelation m_bdry; // Boundary relation from elements to vertices + VertToElemRelation m_cobdry; // Coboundary relation from vertices to elements + + // Mesh Map Definitions + PointMap m_vertexPositions; // vertex position for each vertex + std::vector m_materialVolumeFractionsElement; // the volume fractions of each material for each element + std::vector m_materialVolumeFractionsVertex; // the volume fractions of each material for each vertex + IntMap m_elementParentIDs; // the ID of the parent element from the original mesh + IntMap m_elementDominantMaterials; // the dominant material of the cell (a processed mesh should have only one material per cell) + IntMap m_shapeTypes; // the int enumerator of what type of shape each element is + + int m_numMaterials; // the number of materials present in the mesh + CellTopologyData m_meshTopology; // the topology/connectivity of the mesh + }; + + //-------------------------------------------------------------------------------- +} +} +#endif diff --git a/src/axom/mir/MIRMeshTypes.hpp b/src/axom/mir/MIRMeshTypes.hpp new file mode 100644 index 0000000000..7ec2c79ba2 --- /dev/null +++ b/src/axom/mir/MIRMeshTypes.hpp @@ -0,0 +1,64 @@ +// 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) + +/** + * \file MIRMeshTypes.hpp + * + * \brief Contains the specifications for types aliases used throughout the MIR component. + */ + +#ifndef __MIR_MESH_TYPES_H__ +#define __MIR_MESH_TYPES_H__ + + +#include "axom/core.hpp" // for axom macros +#include "axom/slam.hpp" +#include "axom/primal.hpp" + + +namespace axom +{ +namespace mir +{ + + enum Shape + { + Triangle, + Quad, + Tetrahedron, + Triangular_Prism, + Pyramid, + Hexahedron + }; + + + using Point2 = primal::Point; + + // SET TYPE ALIASES + using PosType = slam::DefaultPositionType; + using ElemType = slam::DefaultElementType; + + using ArrayIndir = slam::policies::ArrayIndirection< PosType, ElemType >; + + using VertSet = slam::PositionSet< PosType, ElemType >; + using ElemSet = slam::PositionSet< PosType, ElemType >; + + // RELATION TYPE ALIASES + using VarCard = slam::policies::VariableCardinality< PosType, ArrayIndir >; + + // Note: This is the actual relation type, which takes in a bunch of policies and data entries to relate to each other. + // Note: It is the relation of the elements to the vertices. + + using ElemToVertRelation = slam::StaticRelation< PosType, ElemType, VarCard, ArrayIndir, ElemSet, VertSet >; + using VertToElemRelation = slam::StaticRelation< PosType, ElemType, VarCard, ArrayIndir, VertSet, ElemSet >; + + // MAP TYPE ALIASES + using BaseSet = slam::Set< PosType, ElemType >; + using ScalarMap = slam::Map< BaseSet, axom::float64 >; + using PointMap = slam::Map< BaseSet, Point2 >; + using IntMap = slam::Map< BaseSet, int >; +} +} +#endif diff --git a/src/axom/mir/MIRUtilities.hpp b/src/axom/mir/MIRUtilities.hpp new file mode 100644 index 0000000000..69c5eb3c17 --- /dev/null +++ b/src/axom/mir/MIRUtilities.hpp @@ -0,0 +1,231 @@ +// 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) + +/** + * \file MIRUtilities.hpp + * + * \brief Contains a set of functions that provide general utility + * within Axom's MIR component. + * + */ + +#ifndef __MIR_UTILITIES_H__ +#define __MIR_UTILITIES_H__ + +#include "ZooClippingTables.hpp" + +//-------------------------------------------------------------------------------- + +namespace axom +{ +namespace mir +{ +namespace utilities +{ + +//-------------------------------------------------------------------------------- + + /** + * \brief Determines the number of vertices of the given shape in the finite element zoo. + * + * \param shape The shape type from the finite element zoo. + * + * \return THe number of vertices of the shape. + */ + inline int numVerts(mir::Shape shape) + { + int numVertices = -1; + switch (shape) + { + case mir::Shape::Triangle: + numVertices = 3; + break; + case mir::Shape::Quad: + numVertices = 4; + break; + case mir::Shape::Tetrahedron: + numVertices = 4; + break; + case mir::Shape::Pyramid: + numVertices = 5; + break; + case mir::Shape::Triangular_Prism: + numVertices = 6; + break; + case mir::Shape::Hexahedron: + numVertices = 8; + break; + default: + printf("Invalid shape. Cannot determine numVerts().\n"); + } + return numVertices; + } + +//-------------------------------------------------------------------------------- + + /** + * \brief Determines the maximum number of possible vertices of the given shape in the finite element zoo. + * This number includes the midpoint vertices between each of the original shape's vertices. + * + * \param shape The shape type from the finite element zoo. + * + * \return THe number of vertices of the shape. + */ + inline int maxPossibleNumVerts(mir::Shape shape) + { + int numVertices = -1; + switch (shape) + { + case mir::Shape::Triangle: + numVertices = 6; + break; + case mir::Shape::Quad: + numVertices = 8; + break; + case mir::Shape::Tetrahedron: + numVertices = 10; + break; + case mir::Shape::Pyramid: + numVertices = 13; + break; + case mir::Shape::Triangular_Prism: + numVertices = 15; + break; + case mir::Shape::Hexahedron: + numVertices = 20; + break; + default: + printf("Invalid shape. Cannot determine maxPossibleNumVerts().\n"); + } + return numVertices; + } + + +//-------------------------------------------------------------------------------- + + /** + * \brief Returns the local vertex ID of the from/to vertex that is one of + * the two endpoints that the edge the given midpoint is on. + * + * \param shapeType The shape type from the finite element zoo. + * \param midpointVertexID The ID of the vertex between the two endpoints. + * \param isFromVertex A flag denoting which of the two edge endpoints to return. + * + * \return The vertex ID of one of the endpoints. + */ + inline int getEdgeEndpoint(const mir::Shape shapeType, + const int midpointVertexID, + const bool isFromVertex) + { + switch(shapeType) + { + case mir::Shape::Triangle: + if ( midpointVertexID == 3 && isFromVertex ) { return 0; } + if ( midpointVertexID == 3 && !isFromVertex ) { return 1; } + if ( midpointVertexID == 4 && isFromVertex ) { return 1; } + if ( midpointVertexID == 4 && !isFromVertex ) { return 2; } + if ( midpointVertexID == 5 && isFromVertex ) { return 2; } + if ( midpointVertexID == 5 && !isFromVertex ) { return 0; } + break; + case mir::Shape::Quad: + if ( midpointVertexID == 4 && isFromVertex ) { return 0; } + if ( midpointVertexID == 4 && !isFromVertex ) { return 1; } + if ( midpointVertexID == 5 && isFromVertex ) { return 1; } + if ( midpointVertexID == 5 && !isFromVertex ) { return 2; } + if ( midpointVertexID == 6 && isFromVertex ) { return 2; } + if ( midpointVertexID == 6 && !isFromVertex ) { return 3; } + if ( midpointVertexID == 7 && isFromVertex ) { return 3; } + if ( midpointVertexID == 7 && !isFromVertex ) { return 0; } + break; + case mir::Shape::Tetrahedron: + if ( midpointVertexID == 4 && isFromVertex ) { return 0; } + if ( midpointVertexID == 4 && !isFromVertex ) { return 1; } + if ( midpointVertexID == 5 && isFromVertex ) { return 0; } + if ( midpointVertexID == 5 && !isFromVertex ) { return 2; } + if ( midpointVertexID == 6 && isFromVertex ) { return 0; } + if ( midpointVertexID == 6 && !isFromVertex ) { return 3; } + if ( midpointVertexID == 7 && isFromVertex ) { return 1; } + if ( midpointVertexID == 7 && !isFromVertex ) { return 2; } + if ( midpointVertexID == 8 && isFromVertex ) { return 2; } + if ( midpointVertexID == 8 && !isFromVertex ) { return 3; } + if ( midpointVertexID == 9 && isFromVertex ) { return 3; } + if ( midpointVertexID == 9 && !isFromVertex ) { return 1; } + case mir::Shape::Pyramid: + if ( midpointVertexID == 5 && isFromVertex ) { return 0; } + if ( midpointVertexID == 5 && !isFromVertex ) { return 1; } + if ( midpointVertexID == 6 && isFromVertex ) { return 0; } + if ( midpointVertexID == 6 && !isFromVertex ) { return 2; } + if ( midpointVertexID == 7 && isFromVertex ) { return 0; } + if ( midpointVertexID == 7 && !isFromVertex ) { return 3; } + if ( midpointVertexID == 8 && isFromVertex ) { return 0; } + if ( midpointVertexID == 8 && !isFromVertex ) { return 4; } + if ( midpointVertexID == 9 && isFromVertex ) { return 1; } + if ( midpointVertexID == 9 && !isFromVertex ) { return 2; } + if ( midpointVertexID == 10 && isFromVertex ) { return 2; } + if ( midpointVertexID == 10 && !isFromVertex ) { return 3; } + if ( midpointVertexID == 11 && isFromVertex ) { return 3; } + if ( midpointVertexID == 11 && !isFromVertex ) { return 4; } + if ( midpointVertexID == 12 && isFromVertex ) { return 4; } + if ( midpointVertexID == 12 && !isFromVertex ) { return 1; } + case mir::Shape::Triangular_Prism: + if ( midpointVertexID == 6 && isFromVertex ) { return 0; } + if ( midpointVertexID == 6 && !isFromVertex ) { return 1; } + if ( midpointVertexID == 7 && isFromVertex ) { return 1; } + if ( midpointVertexID == 7 && !isFromVertex ) { return 2; } + if ( midpointVertexID == 8 && isFromVertex ) { return 2; } + if ( midpointVertexID == 8 && !isFromVertex ) { return 0; } + if ( midpointVertexID == 9 && isFromVertex ) { return 0; } + if ( midpointVertexID == 9 && !isFromVertex ) { return 3; } + if ( midpointVertexID == 10 && isFromVertex ) { return 1; } + if ( midpointVertexID == 10 && !isFromVertex ) { return 4; } + if ( midpointVertexID == 11 && isFromVertex ) { return 2; } + if ( midpointVertexID == 11 && !isFromVertex ) { return 5; } + if ( midpointVertexID == 12 && isFromVertex ) { return 3; } + if ( midpointVertexID == 12 && !isFromVertex ) { return 4; } + if ( midpointVertexID == 13 && isFromVertex ) { return 4; } + if ( midpointVertexID == 13 && !isFromVertex ) { return 5; } + if ( midpointVertexID == 14 && isFromVertex ) { return 5; } + if ( midpointVertexID == 14 && !isFromVertex ) { return 3; } + case mir::Shape::Hexahedron: + if ( midpointVertexID == 8 && isFromVertex ) { return 0; } + if ( midpointVertexID == 8 && !isFromVertex ) { return 1; } + if ( midpointVertexID == 9 && isFromVertex ) { return 1; } + if ( midpointVertexID == 9 && !isFromVertex ) { return 2; } + if ( midpointVertexID == 10 && isFromVertex ) { return 2; } + if ( midpointVertexID == 10 && !isFromVertex ) { return 3; } + if ( midpointVertexID == 11 && isFromVertex ) { return 3; } + if ( midpointVertexID == 11 && !isFromVertex ) { return 0; } + if ( midpointVertexID == 12 && isFromVertex ) { return 0; } + if ( midpointVertexID == 12 && !isFromVertex ) { return 4; } + if ( midpointVertexID == 13 && isFromVertex ) { return 1; } + if ( midpointVertexID == 13 && !isFromVertex ) { return 5; } + if ( midpointVertexID == 14 && isFromVertex ) { return 2; } + if ( midpointVertexID == 14 && !isFromVertex ) { return 6; } + if ( midpointVertexID == 15 && isFromVertex ) { return 3; } + if ( midpointVertexID == 15 && !isFromVertex ) { return 7; } + if ( midpointVertexID == 16 && isFromVertex ) { return 4; } + if ( midpointVertexID == 16 && !isFromVertex ) { return 5; } + if ( midpointVertexID == 17 && isFromVertex ) { return 5; } + if ( midpointVertexID == 17 && !isFromVertex ) { return 6; } + if ( midpointVertexID == 18 && isFromVertex ) { return 6; } + if ( midpointVertexID == 18 && !isFromVertex ) { return 7; } + if ( midpointVertexID == 19 && isFromVertex ) { return 7; } + if ( midpointVertexID == 19 && !isFromVertex ) { return 4; } + default: + printf("Edge endpoint case not implemented.\n"); + return -1; + break; + } + return -1; + } + + +//-------------------------------------------------------------------------------- + +} +} +} + +#endif diff --git a/src/axom/mir/MeshTester.cpp b/src/axom/mir/MeshTester.cpp new file mode 100644 index 0000000000..2b576c66a5 --- /dev/null +++ b/src/axom/mir/MeshTester.cpp @@ -0,0 +1,801 @@ +// 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 "MeshTester.hpp" + +namespace numerics = axom::numerics; +namespace slam = axom::slam; + +namespace axom +{ +namespace mir +{ + +//-------------------------------------------------------------------------------- + +MIRMesh MeshTester::initTestCaseOne() +{ + mir::CellTopologyData topoData; + mir::CellMapData mapData; + mir::CellData cellData; + VolumeFractions volFracs; + + int numElements = 9; + int numVertices = 16; + mir::VertSet verts(numVertices); // Construct a vertex set with 16 vertices + mir::ElemSet elems(numElements); // Construct an element set with 9 elements + + // Create the mesh connectivity information + topoData.m_evInds = { + 0,4,5,1, // elem 0, card 4, start 0 + 1,5,6,2, // elem 1, card 4, start 4 + 2,6,7,3, // elem 2, card 4, start 8 + 4,8,9,5, // elem 3, card 4, start 12 + 5,9,10,6, // elem 4, card 4, start 16 + 6,10,11,7, // elem 5, card 4, start 20 + 8,12,13,9, // elem 6, card 4, start 24 + 9,13,14,10, // elem 7, card 4, start 28 + 10,14,15,11 // elem 8, card 4, start 32, end 36 + }; + + topoData.m_evBegins = { + 0,4,8,12,16,20,24,28,32,36 + }; + topoData.m_veInds = { + 0, // vert 0, card 1, start 0 + 0,1, // vert 1, card 2, start 1 + 1,2, // vert 2, card 2, start 3 + 2, // vert 3, card 1, start 5 + 0,3, // vert 4, card 2, start 6 + 0,1,3,4, // vert 5, card 4, start 8 + 1,2,4,5, // vert 6, card 4, start 12 + 2,5, // vert 7, card 2, start 16 + 3,6, // vert 8, card 2, start 18 + 3,4,6,7, // vert 9, card 4, start 20 + 4,5,7,8, // vert 10, card 4, start 24 + 5,8, // vert 11, card 2, start 28 + 6, // vert 12, card 1, start 30 + 6,7, // vert 13, card 2, start 31 + 7,8, // vert 14, card 2, start 33 + 8, // vert 15, card 1, start 35, end 36 + }; + topoData.m_veBegins = { + 0,1,3,5,6,8,12,16,18,20,24,28,30,31,33,35,36 + }; + + int numMaterials = 2; + enum { GREEN = 0, BLUE = 1 }; + + volFracs.resize(numMaterials); + + volFracs[GREEN] = {1.0, 1.0, 1.0, 1.0, 0.5, 0.2, 0.2, 0.0, 0.0}; + volFracs[BLUE] = {0.0, 0.0, 0.0, 0.0, 0.5, 0.8, 0.8, 1.0, 1.0}; + + + mapData.m_vertexPositions = + { + mir::Point2::make_point( 0.0, 3.0 ), + mir::Point2::make_point( 1.0, 3.0 ), + mir::Point2::make_point( 2.0, 3.0 ), + mir::Point2::make_point( 3.0, 3.0 ), + + mir::Point2::make_point( 0.0, 2.0 ), + mir::Point2::make_point( 1.0, 2.0 ), + mir::Point2::make_point( 2.0, 2.0 ), + mir::Point2::make_point( 3.0, 2.0 ), + + mir::Point2::make_point( 0.0, 1.0 ), + mir::Point2::make_point( 1.0, 1.0 ), + mir::Point2::make_point( 2.0, 1.0 ), + mir::Point2::make_point( 3.0, 1.0 ), + + mir::Point2::make_point( 0.0, 0.0 ), + mir::Point2::make_point( 1.0, 0.0 ), + mir::Point2::make_point( 2.0, 0.0 ), + mir::Point2::make_point( 3.0, 0.0 ) + }; + + mapData.m_elementDominantMaterials = Vec(numElements, NULL_MAT); + mapData.m_elementParents = { 0,1,2,3,4,5,6,7,8 }; // For the base mesh, the parents are always themselves + mapData.m_shapeTypes = Vec(numElements, mir::Shape::Quad); + + // Build the mesh + mir::MIRMesh testMesh; + testMesh.initializeMesh(verts, elems, numMaterials, topoData, mapData, volFracs); + + return testMesh; +} + +//-------------------------------------------------------------------------------- + +mir::MIRMesh MeshTester::initTestCaseTwo() +{ + mir::CellTopologyData topoData; + mir::CellMapData mapData; + mir::CellData cellData; + VolumeFractions volFracs; + + int numElements = 9; + int numVertices = 16; + mir::VertSet verts(numVertices); // Construct a vertex set with 16 vertices + mir::ElemSet elems(numElements); // Construct an element set with 9 elements + + // Create the mesh connectivity information + topoData.m_evInds = { + 0,4,5,1, // elem 0, card 4, start 0 + 1,5,6,2, // elem 1, card 4, start 4 + 2,6,7,3, // elem 2, card 4, start 8 + 4,8,9,5, // elem 3, card 4, start 12 + 5,9,10,6, // elem 4, card 4, start 16 + 6,10,11,7, // elem 5, card 4, start 20 + 8,12,13,9, // elem 6, card 4, start 24 + 9,13,14,10, // elem 7, card 4, start 28 + 10,14,15,11 // elem 8, card 4, start 32, end 36 + }; + + topoData.m_evBegins = { + 0,4,8,12,16,20,24,28,32,36 + }; + topoData.m_veInds = { + 0, // vert 0, card 1, start 0 + 0,1, // vert 1, card 2, start 1 + 1,2, // vert 2, card 2, start 3 + 2, // vert 3, card 1, start 5 + 0,3, // vert 4, card 2, start 6 + 0,1,3,4, // vert 5, card 4, start 8 + 1,2,4,5, // vert 6, card 4, start 12 + 2,5, // vert 7, card 2, start 16 + 3,6, // vert 8, card 2, start 18 + 3,4,6,7, // vert 9, card 4, start 20 + 4,5,7,8, // vert 10, card 4, start 24 + 5,8, // vert 11, card 2, start 28 + 6, // vert 12, card 1, start 30 + 6,7, // vert 13, card 2, start 31 + 7,8, // vert 14, card 2, start 33 + 8, // vert 15, card 1, start 35, end 36 + }; + topoData.m_veBegins = { + 0,1,3,5,6,8,12,16,18,20,24,28,30,31,33,35,36 + }; + + int numMaterials = 3; + enum { BLUE = 0, RED = 1, ORANGE = 2 }; + + volFracs.resize(numMaterials); + volFracs[BLUE] = {1.0, 1.0, 1.0, 1.0, 0.5, 0.2, 0.2, 0.0, 0.0}; + volFracs[RED] = {0.0, 0.0, 0.0, 0.0, 0.3, 0.8, 0.0, 0.3, 1.0}; + volFracs[ORANGE] = {0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.8, 0.7, 0.0}; + + mapData.m_vertexPositions = + { + mir::Point2::make_point( 0.0, 3.0 ), + mir::Point2::make_point( 1.0, 3.0 ), + mir::Point2::make_point( 2.0, 3.0 ), + mir::Point2::make_point( 3.0, 3.0 ), + + mir::Point2::make_point( 0.0, 2.0 ), + mir::Point2::make_point( 1.0, 2.0 ), + mir::Point2::make_point( 2.0, 2.0 ), + mir::Point2::make_point( 3.0, 2.0 ), + + mir::Point2::make_point( 0.0, 1.0 ), + mir::Point2::make_point( 1.0, 1.0 ), + mir::Point2::make_point( 2.0, 1.0 ), + mir::Point2::make_point( 3.0, 1.0 ), + + mir::Point2::make_point( 0.0, 0.0 ), + mir::Point2::make_point( 1.0, 0.0 ), + mir::Point2::make_point( 2.0, 0.0 ), + mir::Point2::make_point( 3.0, 0.0 ) + }; + + mapData.m_elementDominantMaterials = Vec(numElements, NULL_MAT); + mapData.m_elementParents = { 0,1,2,3,4,5,6,7,8 }; // For the base mesh, the parents are always themselves + mapData.m_shapeTypes = Vec(numElements, mir::Shape::Quad); + + // Build the mesh + mir::MIRMesh testMesh; + testMesh.initializeMesh(verts, elems, numMaterials, topoData, mapData, volFracs); + + return testMesh; +} + +//-------------------------------------------------------------------------------- + +mir::MIRMesh MeshTester::initTestCaseThree() +{ + mir::CellTopologyData topoData; + mir::CellMapData mapData; + mir::CellData cellData; + VolumeFractions volFracs; + + int numElements = 4; + int numVertices = 6; // OR create a middle triangle with all of one material, and then a ring of triangles around it that are full of the other material + + mir::VertSet verts = mir::VertSet(numVertices); + mir::ElemSet elems = mir::ElemSet(numElements); + + // Create the mesh connectivity information + topoData.m_evInds = { + 0,1,2, // elem 0, card 3, start 0 + 1,3,4, // elem 1, card 3, start 3 + 1,4,2, // elem 2, card 3, start 6 + 2,4,5 // elem 3, card 3, start 9, end 12 + }; + + topoData.m_evBegins = { + 0,3,6,9,12 + }; + topoData.m_veInds = { + 0, // vert 0, card 1, start 0 + 0,1,2, // vert 1, card 3, start 1 + 0,2,3, // vert 2, card 3, start 4 + 1, // vert 3, card 1, start 7 + 1,2,3, // vert 4, card 3, start 8 + 3 // vert 5, card 1, start 11, end 12 + }; + topoData.m_veBegins = { + 0,1,4,7,8,11,12 + }; + + + int numMaterials = 2; + enum { BLUE = 0, RED = 1, }; + + volFracs.resize(numMaterials); + volFracs[BLUE] = {0.0, 0.5, 0.8, 0.5}; + volFracs[RED] = {1.0, 0.5, 0.2, 0.5}; + + mapData.m_vertexPositions = + { + mir::Point2::make_point( 1.0, 2.0 ), + mir::Point2::make_point( 0.5, 1.0 ), + mir::Point2::make_point( 1.5, 1.0 ), + mir::Point2::make_point( 0.0, 0.0 ), + mir::Point2::make_point( 1.0, 0.0 ), + mir::Point2::make_point( 2.0, 0.0 ) + }; + + mapData.m_elementDominantMaterials = Vec(numElements, NULL_MAT); + mapData.m_elementParents = { 0,1,2,3 }; // For the base mesh, the parents are always themselves + mapData.m_shapeTypes = Vec(numElements, mir::Shape::Triangle); + + // Build the mesh + mir::MIRMesh testMesh; + testMesh.initializeMesh(verts, elems, numMaterials, topoData, mapData, volFracs); + + return testMesh; +} + +//-------------------------------------------------------------------------------- + +mir::MIRMesh MeshTester::initTestCaseFour() +{ + mir::CellTopologyData topoData; + mir::CellMapData mapData; + mir::CellData cellData; + VolumeFractions volFracs; + + int numElements = 9; + int numVertices = 16; + mir::VertSet verts = mir::VertSet(numVertices); // Construct a vertex set with 16 vertices + mir::ElemSet elems = mir::ElemSet(numElements); // Construct an element set with 9 elements + + // Create the mesh connectivity information + topoData.m_evInds = { + 0,4,5,1, // elem 0, card 4, start 0 + 1,5,6,2, // elem 1, card 4, start 4 + 2,6,7,3, // elem 2, card 4, start 8 + 4,8,9,5, // elem 3, card 4, start 12 + 5,9,10,6, // elem 4, card 4, start 16 + 6,10,11,7, // elem 5, card 4, start 20 + 8,12,13,9, // elem 6, card 4, start 24 + 9,13,14,10, // elem 7, card 4, start 28 + 10,14,15,11 // elem 8, card 4, start 32, end 36 + }; + + topoData.m_evBegins = { + 0,4,8,12,16,20,24,28,32,36 + }; + topoData.m_veInds = { + 0, // vert 0, card 1, start 0 + 0,1, // vert 1, card 2, start 1 + 1,2, // vert 2, card 2, start 3 + 2, // vert 3, card 1, start 5 + 0,3, // vert 4, card 2, start 6 + 0,1,3,4, // vert 5, card 4, start 8 + 1,2,4,5, // vert 6, card 4, start 12 + 2,5, // vert 7, card 2, start 16 + 3,6, // vert 8, card 2, start 18 + 3,4,6,7, // vert 9, card 4, start 20 + 4,5,7,8, // vert 10, card 4, start 24 + 5,8, // vert 11, card 2, start 28 + 6, // vert 12, card 1, start 30 + 6,7, // vert 13, card 2, start 31 + 7,8, // vert 14, card 2, start 33 + 8, // vert 15, card 1, start 35, end 36 + }; + topoData.m_veBegins = { + 0,1,3,5,6,8,12,16,18,20,24,28,30,31,33,35,36 + }; + + + mapData.m_vertexPositions = + { + mir::Point2::make_point( 0.0, 3.0 ), + mir::Point2::make_point( 1.0, 3.0 ), + mir::Point2::make_point( 2.0, 3.0 ), + mir::Point2::make_point( 3.0, 3.0 ), + + mir::Point2::make_point( 0.0, 2.0 ), + mir::Point2::make_point( 1.0, 2.0 ), + mir::Point2::make_point( 2.0, 2.0 ), + mir::Point2::make_point( 3.0, 2.0 ), + + mir::Point2::make_point( 0.0, 1.0 ), + mir::Point2::make_point( 1.0, 1.0 ), + mir::Point2::make_point( 2.0, 1.0 ), + mir::Point2::make_point( 3.0, 1.0 ), + + mir::Point2::make_point( 0.0, 0.0 ), + mir::Point2::make_point( 1.0, 0.0 ), + mir::Point2::make_point( 2.0, 0.0 ), + mir::Point2::make_point( 3.0, 0.0 ) + }; + + int numMaterials = 2; + enum { GREEN = 0, BLUE = 1 }; + + volFracs.resize(numMaterials); + + auto& greenVolumeFractions = volFracs[GREEN]; + auto& blueVolumeFractions = volFracs[BLUE]; + const auto& points = mapData.m_vertexPositions; + const auto& evInds = topoData.m_evInds; + + greenVolumeFractions.resize(numElements); + blueVolumeFractions.resize(numElements); + + // Generate the element volume fractions for the circle + auto circleCenter = mir::Point2::make_point(1.5, 1.5); + axom::float64 circleRadius = 1.25; + int gridSize = 1000; + for (int i = 0; i < numElements; ++i) + { + auto vf = calculatePercentOverlapMonteCarlo(gridSize, + circleCenter, + circleRadius, + points[evInds[i * 4 + 0]], + points[evInds[i * 4 + 1]], + points[evInds[i * 4 + 2]], + points[evInds[i * 4 + 3]]); + greenVolumeFractions[i] = vf; + blueVolumeFractions[i] = 1.0 - vf; + } + + mapData.m_elementDominantMaterials = Vec(numElements, NULL_MAT); + mapData.m_elementParents = { 0,1,2,3,4,5,6,7,8 }; // For the base mesh, the parents are always themselves + mapData.m_shapeTypes = Vec(numElements, mir::Shape::Quad); + + // Build the mesh + mir::MIRMesh testMesh; + testMesh.initializeMesh(verts, elems, numMaterials, topoData, mapData, volFracs); + + return testMesh; +} + +//-------------------------------------------------------------------------------- + +mir::MIRMesh MeshTester::createUniformGridTestCaseMesh(int gridSize, const mir::Point2& circleCenter, axom::float64 circleRadius) +{ + // Generate the mesh topology + mir::CellData cellData = generateGrid(gridSize); + + mir::VertSet verts = mir::VertSet(cellData.m_numVerts); // Construct the vertex set + mir::ElemSet elems = mir::ElemSet(cellData.m_numElems); // Construct the element set + + int numMaterials = 2; + enum { GREEN = 0, BLUE = 1 }; + + VolumeFractions volFracs; + volFracs.resize(numMaterials); + volFracs[GREEN].resize(cellData.m_numElems); + volFracs[BLUE].resize(cellData.m_numElems); + + // Generate the element volume fractions for the circle + const int numMonteCarloSamples = 100; + auto& pos = cellData.m_mapData.m_vertexPositions; + const auto& evInds = cellData.m_topology.m_evInds; + for (int i = 0; i < cellData.m_numElems; ++i) + { + auto vf = calculatePercentOverlapMonteCarlo(numMonteCarloSamples, + circleCenter, + circleRadius, + pos[evInds[i * 4 + 0]], + pos[evInds[i * 4 + 1]], + pos[evInds[i * 4 + 2]], + pos[evInds[i * 4 + 3]]); + volFracs[GREEN][i] = vf; + volFracs[BLUE][i] = 1.0 - vf; + } + + cellData.m_mapData.m_elementDominantMaterials = Vec(cellData.m_numVerts, NULL_MAT); + cellData.m_mapData.m_shapeTypes = Vec(cellData.m_numVerts, mir::Shape::Quad); + cellData.m_mapData.m_elementParents.resize(cellData.m_numVerts); + for (auto i : elems.positions() ) + { + cellData.m_mapData.m_elementParents[i] = i; + } + + // Build the mesh + mir::MIRMesh testMesh; + testMesh.initializeMesh(verts, elems, numMaterials, cellData.m_topology, cellData.m_mapData, volFracs); + + return testMesh; +} + +//-------------------------------------------------------------------------------- + +axom::float64 MeshTester::calculatePercentOverlapMonteCarlo( + int gridSize, + const mir::Point2& circleCenter, + axom::float64 circleRadius, + const mir::Point2& quadP0, + const mir::Point2& quadP1, + const mir::Point2& quadP2, + const mir::Point2& quadP3) +{ + // Check if any of the quad's corners are within the circle + auto d0Sq = primal::squared_distance(quadP0, circleCenter); + auto d1Sq = primal::squared_distance(quadP1, circleCenter); + auto d2Sq = primal::squared_distance(quadP2, circleCenter); + auto d3Sq = primal::squared_distance(quadP3, circleCenter); + auto dRSq = circleRadius * circleRadius; + + int inFlags = ((d0Sq < dRSq) ? 1 << 0 : 0) + + ((d1Sq < dRSq) ? 1 << 1 : 0) + + ((d2Sq < dRSq) ? 1 << 2 : 0) + + ((d3Sq < dRSq) ? 1 << 3 : 0); + const int allFlags = 15; + const int noFlags = 0; + + if (inFlags == allFlags ) + { + // The entire quad overlaps the circle + return 1.; + } + else if( inFlags == noFlags) + { + return 0.; + } + else + { + // Some of the quad overlaps the circle, so run the Monte Carlo sampling to determine how much + axom::float64 delta_x = axom::utilities::abs(quadP2[0] - quadP1[0]) / static_cast(gridSize - 1); + axom::float64 delta_y = axom::utilities::abs(quadP0[1] - quadP1[1]) / static_cast(gridSize - 1); + int countOverlap = 0; + for (int y = 0; y < gridSize; ++y) + { + for (int x = 0; x < gridSize; ++x) + { + mir::Point2 samplePoint = mir::Point2::make_point(delta_x * x + quadP1[0], + delta_y * y + quadP1[1]); + if (primal::squared_distance(samplePoint, circleCenter) < dRSq) + ++countOverlap; + } + } + return countOverlap / static_cast(gridSize * gridSize); + } +} + +//-------------------------------------------------------------------------------- + +mir::CellData MeshTester::generateGrid(int gridSize) +{ + // Generate the topology for a uniform quad mesh with n x n elements automatically + int numElements = gridSize * gridSize; + int numVertices = (gridSize + 1) * (gridSize + 1); + + mir::CellData data; + + data.m_numVerts = numVertices; + data.m_numElems = numElements; + + // Generate the evInds + auto& evInds = data.m_topology.m_evInds; + for (int eID = 0; eID < numElements; ++eID) + { + int row = eID / gridSize; // note the integer division + int vertsPerRow = gridSize + 1; + int elemsPerRow = gridSize; + + evInds.push_back( (eID % elemsPerRow) + row * vertsPerRow + 0); + evInds.push_back( (eID % elemsPerRow) + (row + 1) * vertsPerRow + 0); + evInds.push_back( (eID % elemsPerRow) + (row + 1) * vertsPerRow + 1); + evInds.push_back( (eID % elemsPerRow) + row * vertsPerRow + 1); + } + + // Generate the evBegins + auto& evBegins = data.m_topology.m_evBegins; + evBegins.push_back(0); + for (int i = 0; i < numElements; ++i) + { + evBegins.push_back((i + 1) * 4); + } + + // Generate the veInds + auto& veInds = data.m_topology.m_veInds; + auto& veBegins = data.m_topology.m_veBegins; + std::map > veInds_data; + for (int evInd_itr = 0; evInd_itr < numElements * 4; ++evInd_itr) + { + int currentElementID = evInd_itr / 4; // note the integer division + veInds_data[evInds[evInd_itr]].push_back(currentElementID); + } + + for (auto itr = veInds_data.begin(); itr != veInds_data.end(); itr++) + { + // Sort the vector + std::sort(itr->second.begin(), itr->second.end()); + + // Add the elements associated with the current vertex to veInds + for (unsigned long i = 0; i < itr->second.size(); ++i) + veInds.push_back(itr->second[i]); + } + + // Generate the veBegins + veBegins.push_back(0); + int currentIndexCount = 0; + for (auto itr = veInds_data.begin(); itr != veInds_data.end(); itr++) + { + currentIndexCount += itr->second.size(); + veBegins.push_back(currentIndexCount); + } + + // Generate the vertex positions + auto& points = data.m_mapData.m_vertexPositions; + for (int y = gridSize; y > -1; --y) + { + for (int x = 0; x < gridSize + 1; ++x) + { + points.push_back(mir::Point2::make_point(x, y)); + } + } + + // // Print out the results + // printf("evInds: { "); + // for (int i = 0; i < evInds.size(); i++) + // { + // printf("%d ", evInds[i]); + // if ((i+1) % 4 == 0 && i != 0) + // printf("\n"); + // } + // printf("}\n"); + + // printf("evBegins: { "); + // for (int i = 0; i < evBegins.size(); i++) + // { + // printf("%d ", evBegins[i]); + // } + // printf("}\n"); + + // printf("veInds: { "); + // for (int i = 0; i < veInds.size(); i++) + // { + // printf("%d ", veInds[i]); + // } + // printf("}\n"); + + // printf("veBegins: { "); + // for (int i = 0; i < veBegins.size(); i++) + // { + // printf("%d ", veBegins[i]); + // } + // printf("}\n"); + + // printf("points: { "); + // for (int i = 0; i < numVertices; ++i) + // { + // printf("{%.2f, %.2f} ", points[i][0], points[i][1]); + // } + // printf("}\n"); + + return data; +} + +//-------------------------------------------------------------------------------- + +mir::MIRMesh MeshTester::initTestCaseFive(int gridSize, int numCircles) +{ + + // Generate the mesh topology + mir::CellData cellData = generateGrid(gridSize); + + mir::VertSet verts = mir::VertSet(cellData.m_numVerts); // Construct the vertex set + mir::ElemSet elems = mir::ElemSet(cellData.m_numElems); // Construct the element set + + // Generate the element volume fractions with concentric circles + int numMaterials = numCircles + 1; + int defaultMaterialID = numMaterials - 1; // default material is always the last index + + mir::Point2 circleCenter = mir::Point2::make_point(gridSize / 2.0, gridSize / 2.0); // all circles are centered around the same point + + // Initialize the radii of the circles + std::vector circleRadii; + axom::float64 maxRadius = gridSize / 2.4; // Note: The choice of divisor is arbitrary + axom::float64 minRadius = gridSize / 8; // Note: The choice of divisor is arbitrary + + axom::float64 radiusDelta; + if (numCircles <= 1) + radiusDelta = (maxRadius - minRadius); + else + radiusDelta = (maxRadius - minRadius) / (double) (numCircles - 1); + + for (int i = 0; i < numCircles; ++i) + { + circleRadii.push_back( minRadius + (i * radiusDelta) ); + } + + // Initialize all material volume fractions to 0 + std::vector > materialVolumeFractionsData; + for (int i = 0; i < numMaterials; ++i) + { + std::vector tempVec; + tempVec.resize(cellData.m_numElems); + materialVolumeFractionsData.push_back(tempVec); + } + + // Use the uniform sampling method to generate volume fractions for each material + // Note: Assumes that the cell is a parallelogram. This could be modified via biliear interpolation + for (int eID = 0; eID < cellData.m_numElems; ++eID) + { + mir::Point2& v0 = cellData.m_mapData.m_vertexPositions[cellData.m_topology.m_evInds[eID * 4 + 0]]; + mir::Point2& v1 = cellData.m_mapData.m_vertexPositions[cellData.m_topology.m_evInds[eID * 4 + 1]]; + mir::Point2& v2 = cellData.m_mapData.m_vertexPositions[cellData.m_topology.m_evInds[eID * 4 + 2]]; + //mir::Point2& v3 = cellData.m_mapData.m_vertexPositions[cellData.m_topology.m_evInds[eID * 4 + 3]]; + + // Run the uniform sampling to determine how much of the current cell is composed of each material + int materialCount[numMaterials]; for (int i = 0; i < numMaterials; ++i) materialCount[i] = 0; + + for (int matID = 0; matID < numMaterials; ++matID) + { + materialVolumeFractionsData[matID][eID] = materialCount[matID] / (double) (gridSize * gridSize); + } + + axom::float64 delta_x = axom::utilities::abs(v2[0] - v1[0]) / (double) (gridSize - 1); + axom::float64 delta_y = axom::utilities::abs(v0[1] - v1[1]) / (double) (gridSize - 1); + + for (int y = 0; y < gridSize; ++y) + { + for (int x = 0; x < gridSize; ++x) + { + mir::Point2 samplePoint = mir::Point2::make_point(delta_x * x + v1[0], delta_y * y + v1[1]); + bool isPointSampled = false; + for (int cID = 0; cID < numCircles && !isPointSampled; ++cID) + { + const auto r = circleRadii[cID]; + if (primal::squared_distance(samplePoint, circleCenter) < r*r) + { + materialCount[cID]++; + isPointSampled = true; + } + } + if (!isPointSampled) + { + // The point was not within any of the circles, so increment the count for the default material + materialCount[defaultMaterialID]++; + } + } + } + + // Assign the element volume fractions based on the count of the samples in each circle + for (int matID = 0; matID < numMaterials; ++matID) + { + materialVolumeFractionsData[matID][eID] = materialCount[matID] / (double) (gridSize * gridSize); + } + } + + std::vector elementParents; // For the base mesh, the parents are always themselves + std::vector elementDominantMaterials; + std::vector elementShapeTypes; + for (int i = 0; i < cellData.m_numElems; ++i) + { + elementParents.push_back(i); + elementDominantMaterials.push_back(NULL_MAT); + elementShapeTypes.push_back(mir::Shape::Quad); + } + + CellTopologyData topology; + topology.m_evInds = cellData.m_topology.m_evInds; + topology.m_evBegins = cellData.m_topology.m_evBegins; + topology.m_veInds = cellData.m_topology.m_veInds; + topology.m_veBegins = cellData.m_topology.m_veBegins; + + CellMapData mapData; + mapData.m_elementDominantMaterials = elementDominantMaterials; + mapData.m_elementParents = elementParents; + mapData.m_vertexPositions = cellData.m_mapData.m_vertexPositions; + mapData.m_shapeTypes = elementShapeTypes; + + // Build the mesh + mir::MIRMesh testMesh; + testMesh.initializeMesh(verts, elems, numMaterials, topology, mapData, materialVolumeFractionsData); + + return testMesh; +} + +//-------------------------------------------------------------------------------- + +int MeshTester::circleQuadCornersOverlaps(const mir::Point2& circleCenter, + axom::float64 circleRadius, + const mir::Point2& quadP0, + const mir::Point2& quadP1, + const mir::Point2& quadP2, + const mir::Point2& quadP3) +{ + // Check if any of the quad's corners are within the circle + auto d0Sq = primal::squared_distance(quadP0, circleCenter); + auto d1Sq = primal::squared_distance(quadP1, circleCenter); + auto d2Sq = primal::squared_distance(quadP2, circleCenter); + auto d3Sq = primal::squared_distance(quadP3, circleCenter); + auto dRSq = circleRadius * circleRadius; + + int numCorners = 0; + + if (d0Sq < dRSq) + numCorners++; + if (d1Sq < dRSq) + numCorners++; + if (d2Sq < dRSq) + numCorners++; + if (d3Sq < dRSq) + numCorners++; + + return numCorners; +} + +//-------------------------------------------------------------------------------- + +mir::MIRMesh MeshTester::initQuadClippingTestMesh() +{ + // Generate the mesh topology + int gridSize = 3; + mir::CellData cellData = generateGrid(gridSize); + + mir::VertSet verts = mir::VertSet(cellData.m_numVerts); // Construct the vertex set + mir::ElemSet elems = mir::ElemSet(cellData.m_numElems); // Construct the element set + + int numMaterials = 2; + + std::vector > elementVF; + elementVF.resize(numMaterials); + elementVF[0] = {1.0, 1.0, 1.0, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0}; + elementVF[1] = {0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0}; + + std::vector elementParents; + std::vector elementDominantMaterials; + std::vector elementShapeTypes; + for (int i = 0; i < cellData.m_numElems; ++i) + { + elementParents.push_back(i); + elementDominantMaterials.push_back(NULL_MAT); + elementShapeTypes.push_back(mir::Shape::Quad); + } + + cellData.m_mapData.m_elementDominantMaterials = elementDominantMaterials; + cellData.m_mapData.m_elementParents = elementParents; + cellData.m_mapData.m_shapeTypes = elementShapeTypes; + + // Build the mesh + mir::MIRMesh testMesh; + testMesh.initializeMesh(verts, elems, numMaterials, cellData.m_topology, cellData.m_mapData, elementVF); + + return testMesh; +} + +//-------------------------------------------------------------------------------- + +} +} diff --git a/src/axom/mir/MeshTester.hpp b/src/axom/mir/MeshTester.hpp new file mode 100644 index 0000000000..dcbc380460 --- /dev/null +++ b/src/axom/mir/MeshTester.hpp @@ -0,0 +1,184 @@ +// 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) + +/** + * \file MeshTester.hpp + * + * \brief Contains the specification for the MeshTester class. + * + */ + +#ifndef __MESH_TESTER_H__ +#define __MESH_TESTER_H__ + +#include "axom/core.hpp" // for axom macros +#include "axom/slam.hpp" // unified header for slam classes and functions + +#include "MIRMesh.hpp" +#include "MIRUtilities.hpp" + +#include +#include + +namespace axom +{ +namespace mir +{ + /** + * \class MeshTester + * + * \brief A class used to generate MIRMeshs with specific properties so + * that the reconstruction output can be validated visually. + * + */ + class MeshTester + { +public: + template + using Vec = std::vector; + + using IndexVec = Vec; + using VolFracVec = Vec; + using VolumeFractions = Vec; + + + public: + /** + * \brief Default constructor. + */ + MeshTester() = default; + + /** + * \brief Default destructor. + */ + ~MeshTester() = default; + + public: + /** + * \brief Initializes an MIRMesh based on the example from Meredith 2004 paper. + * + * \note The mesh is a 3x3 uniform grid of quads with 2 materials. + * + * \return The generated mesh. + */ + mir::MIRMesh initTestCaseOne(); + + /** + * \brief Initializes an MIRMesh based on the example from Meredith and Childs 2010 paper. + * + * \note The mesh is a 3x3 uniform grid of quads with 3 materials. + * + * \return The generated mesh. + */ + mir::MIRMesh initTestCaseTwo(); + + /** + * \brief Initializes an MIRMesh used for testing triangle clipping cases. + * + * \note The mesh is a set of four triangles with 2 materials + * + * \return The generated mesh. + */ + mir::MIRMesh initTestCaseThree(); + + /** + * \brief Intializes a mesh used for testing a single circle of one materials surrounded by another. + * + * \note The mesh is a 3x3 uniform grid with 2 materials and has a single circle in the center composed + * of one material and is surrounded by a second material. + * + * \return The generated mesh. + */ + mir::MIRMesh initTestCaseFour(); + + /** + * \brief Initializes a mesh to be used for testing a set of concentric circles centered in a uniform grid. + * + * \param gridSize The number of elements in the width and the height of the uniform grid. + * \param numCircles The number of concentric circles that are centered in the grid. + * + * \note Each circle is composed of a different material. + * + * \return The generated mesh. + */ + mir::MIRMesh initTestCaseFive(int gridSize, int numCircles); // multiple materials, multiple concentric circles + + /** + * \brief Initializes a mesh composed of a uniform grid with a circle of material in it. + * + * \param gridSize The number of elements in the width and height of the uniform grid. + * \param circleCenter The center point of the circle. + * \param circleRadius The radius of the circle. + * + * \return The generated mesh. + */ + mir::MIRMesh createUniformGridTestCaseMesh(int gridSize, + const mir::Point2& circleCenter, + axom::float64 circleRadius); + + /** + * \brief Initializes a mesh to be used for validating the results of quad clipping. + * + * \note The mesh is a 3x3 uniform grid with 2 materials and element volume fraction such + * that the mesh would be split horizontally through the middle. + * + * \return The generated mesh. + */ + mir::MIRMesh initQuadClippingTestMesh(); + + private: + /** + * \brief Generates a 2D uniform grid of n x n elements. + * + * \param gridSize The number of elements in the width and height of the uniform grid. + */ + mir::CellData generateGrid(int gridSize); + + /** + * \brief Calculates the percent overlap between the given circle and quad. + * + * \param gridSize The size of the uniform grid which will be sampled over to check for overlap. + * \param circleCenter The center point of the circle. + * \param circleRadius The radius of the circle. + * \param quadP0 The upper left vertex of the quad. + * \param quadP1 The lower left vertex of the quad. + * \param quadP2 The lower right vertex of the quad. + * \param quadP3 The upper right vertex of the quad. + * + * /return The percent value overlap of the circle and the quad between [0, 1]. + */ + axom::float64 calculatePercentOverlapMonteCarlo( + int gridSize, + const mir::Point2& circleCenter, + axom::float64 circleRadius, + const mir::Point2& quadP0, + const mir::Point2& quadP1, + const mir::Point2& quadP2, + const mir::Point2& quadP3); + + /** + * \brief Calculates the number of corners of the quad that are within the circle. + * + * \param circleCenter The center point of the circle. + * \param circleRadius The radius of the circle. + * \param quadP0 The upper left vertex of the quad. + * \param quadP1 The lower left vertex of the quad. + * \param quadP2 The lower right vertex of the quad. + * \param quadP3 The upper right vertex of the quad. + * + * \return The number of corners of the quad that are within the circle. + */ + int circleQuadCornersOverlaps(const mir::Point2& circleCenter, + axom::float64 circleRadius, + const mir::Point2& quadP0, + const mir::Point2& quadP1, + const mir::Point2& quadP2, + const mir::Point2& quadP3); + + }; +} +} + +#endif diff --git a/src/axom/mir/README.md b/src/axom/mir/README.md new file mode 100644 index 0000000000..fa6b7fba6a --- /dev/null +++ b/src/axom/mir/README.md @@ -0,0 +1,8 @@ +MIR: Material Interface Reconstruction {#mirtop} +================================================ + +[MIR](@ref axom::mir) provides algorithms for reconstructing interfaces +between materials in a multimaterial mesh. + +The [Mir user documentation](../../../sphinx/axom_docs/html/axom/mir/docs/sphinx/index.html) +describes these algorithms. diff --git a/src/axom/mir/ZooClippingTables.cpp b/src/axom/mir/ZooClippingTables.cpp new file mode 100644 index 0000000000..9f87edbe5d --- /dev/null +++ b/src/axom/mir/ZooClippingTables.cpp @@ -0,0 +1,104 @@ +// 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 "ZooClippingTables.hpp" + +namespace axom +{ +namespace mir +{ + + // Quad Vertex Local Indices + // + // 0 7 3 + // @---------@---------@ + // | | + // | | + // | | + // 4 @ @ 6 + // | | + // | | + // | | + // @---------@---------@ + // 1 5 2 + // + const int quadClipTable[16][19] = + { + {4,0,1,2,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, + {3,3,7,6,4,7,0,2,6,3,0,1,2,-1,-1,-1,-1,-1,-1}, + {3,6,5,2,4,3,1,5,6,3,0,1,3,-1,-1,-1,-1,-1,-1}, + {4,0,1,5,7,4,7,5,2,3,-1,-1,-1,-1,-1,-1,-1,-1}, + {3,4,1,5,4,0,4,5,2,3,0,2,3,-1,-1,-1,-1,-1,-1}, + {3,4,1,5,4,0,4,5,2,4,0,2,6,7,3,7,6,3,-1}, + {4,0,4,6,3,4,4,1,2,6,-1,-1,-1,-1,-1,-1,-1,-1,-1}, + {3,0,4,7,4,4,1,3,7,3,1,2,3,-1,-1,-1,-1,-1,-1}, + {3,0,4,7,4,4,1,3,7,3,1,2,3,-1,-1,-1,-1,-1,-1}, + {4,0,4,6,3,4,4,1,2,6,-1,-1,-1,-1,-1,-1,-1,-1,-1}, + {3,4,1,5,4,0,4,5,2,4,0,2,6,7,3,7,6,3,-1}, + {3,4,1,5,4,0,4,5,2,3,0,2,3,-1,-1,-1,-1,-1,-1}, + {4,0,1,5,7,4,7,5,2,3,-1,-1,-1,-1,-1,-1,-1,-1}, + {3,6,5,2,4,3,1,5,6,3,0,1,3,-1,-1,-1,-1,-1,-1}, + {3,3,7,6,4,7,0,2,6,3,0,1,2,-1,-1,-1,-1,-1,-1}, + {4,0,1,2,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1} + }; + + const std::vector > quadClipTableVec = + { + {4,0,1,2,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, + {3,3,7,6,4,7,0,2,6,3,0,1,2,-1,-1,-1,-1,-1,-1}, + {3,6,5,2,4,3,1,5,6,3,0,1,3,-1,-1,-1,-1,-1,-1}, + {4,0,1,5,7,4,7,5,2,3,-1,-1,-1,-1,-1,-1,-1,-1}, + {3,4,1,5,4,0,4,5,2,3,0,2,3,-1,-1,-1,-1,-1,-1}, + {3,4,1,5,4,0,4,5,2,4,0,2,6,7,3,7,6,3,-1}, + {4,0,4,6,3,4,4,1,2,6,-1,-1,-1,-1,-1,-1,-1,-1,-1}, + {3,0,4,7,4,4,1,3,7,3,1,2,3,-1,-1,-1,-1,-1,-1}, + {3,0,4,7,4,4,1,3,7,3,1,2,3,-1,-1,-1,-1,-1,-1}, + {4,0,4,6,3,4,4,1,2,6,-1,-1,-1,-1,-1,-1,-1,-1,-1}, + {3,4,1,5,4,0,4,5,2,4,0,2,6,7,3,7,6,3,-1}, + {3,4,1,5,4,0,4,5,2,3,0,2,3,-1,-1,-1,-1,-1,-1}, + {4,0,1,5,7,4,7,5,2,3,-1,-1,-1,-1,-1,-1,-1,-1}, + {3,6,5,2,4,3,1,5,6,3,0,1,3,-1,-1,-1,-1,-1,-1}, + {3,3,7,6,4,7,0,2,6,3,0,1,2,-1,-1,-1,-1,-1,-1}, + {4,0,1,2,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1} + }; + + + + // Triangle Vertex Local Indices + // 0 + // @ + // / \ + // / \ + // 3 @ @ 5 + // / \ + // / \ + // 1 @-----@-----@ 2 + // 4 + const int triangleClipTable[8][10] = + { + {3,0,1,2,-1,-1,-1,-1,-1,-1}, + {4,0,1,4,5,3,5,4,2,-1}, + {4,0,3,4,2,3,3,1,4,-1}, + {4,3,1,2,5,3,0,3,5,-1}, + {4,3,1,2,5,3,0,3,5,-1}, + {4,0,3,4,2,3,3,1,4,-1}, + {4,0,1,4,5,3,5,4,2,-1}, + {3,0,1,2,-1,-1,-1,-1,-1,-1} + }; + + const std::vector > triangleClipTableVec = + { + {3,0,1,2,-1,-1,-1,-1,-1,-1}, + {4,0,1,4,5,3,5,4,2,-1}, + {4,0,3,4,2,3,3,1,4,-1}, + {4,3,1,2,5,3,0,3,5,-1}, + {4,3,1,2,5,3,0,3,5,-1}, + {4,0,3,4,2,3,3,1,4,-1}, + {4,0,1,4,5,3,5,4,2,-1}, + {3,0,1,2,-1,-1,-1,-1,-1,-1} + }; + +} +} \ No newline at end of file diff --git a/src/axom/mir/ZooClippingTables.hpp b/src/axom/mir/ZooClippingTables.hpp new file mode 100644 index 0000000000..de2b8be83c --- /dev/null +++ b/src/axom/mir/ZooClippingTables.hpp @@ -0,0 +1,29 @@ +// 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) + +#ifndef __ZOO_CLIPPING_TABLES_H__ +#define __ZOO_CLIPPING_TABLES_H__ + +/** + * \file ZooClippingTables.hpp + * + * \brief Contains the definitions for the clipping cases and enumerator + * for the shape types in the zoo. + */ + +#include + +namespace axom +{ +namespace mir +{ + extern const int quadClipTable[16][19]; + extern const int triangleClipTable[8][10]; + + extern const std::vector > triangleClipTableVec; + extern const std::vector > quadClipTableVec; +} +} +#endif diff --git a/src/axom/mir/docs/sphinx/equi_z.rst b/src/axom/mir/docs/sphinx/equi_z.rst new file mode 100644 index 0000000000..3201fb6660 --- /dev/null +++ b/src/axom/mir/docs/sphinx/equi_z.rst @@ -0,0 +1,61 @@ +.. ## 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) + +======================= +Equi-Z Algorithm +======================= + +The Equi-Z algorithm provided here is an implementation of the algorithm described in +the 2005 paper "Material Interface Reconstruction in VisIt" by J. S. Meredith. + +Functionality for this algorithm is provided by ``computeReconstructedInterface()``. +It takes in references to two meshes: the mesh to be processed and the mesh that +will contain the results of the material interface reconstruction. + +.. raw:: html + +

Key Concepts

+ +Full details of the algorithm can be found in the original paper, but key concepts +are explained here. + +The input into this algorithm is a volumetric mesh composed of element based volume fractions. +The elements of this input mesh can contain multiple materials. The output of this algorithm +is a new volumetric mesh composed of clean elements (i.e. only one material is present in each element). +Note that both the input and output meshes are composed of elements that are part of the finite element zoo (described below). +Each element of the mesh is then processed independently. For each pair of materials of the current element, convert the element +centered volume fractions to vertex centered volume fractions, which are found by calculating the average value of the materials +of the adjacent elements. Then, determine the dominant material at each vertex of the element. The configuration of the dominant +materials is used as an index into a lookup table of clipping cases (described below). Once the clipping case is found, linear interpolation is used +to find the exact points along the edges at which the current element should be split. The element is then decomposed into clean, zoo +elements. + + +.. raw:: html + +

Finite Element Zoo

+ +For two-dimensional meshes, the finite element zoo consists of triangles and quadrilaterals, while for three-dimensional meshes, +the finite element zoo consists of tetrahedrons, pyramids, triangular prisms (wedges), and hexahedrons. This can be seen in the +figure below. + +.. figure:: zoo.png + +Splitting the mesh into elements contained in the zoo prevents a couple of issues: the first being +interpolation imprinting and the second being the need to split arbitrary polyhedra. + +.. raw:: html + +

Clipping Cases

+ +During the course of the algorithm, it is necessary to map each shape and the material that is dominant at each of its vertices +to a specific, pre-defined clipping case. For two dimensional shapes, the cases can be seen below. + +.. figure:: triangle_cases.png + +.. figure:: quad_cases.png + +Note that while this algorithm can handle any number of materials in the input mesh, it considers two materials at a time when splitting +an element. \ No newline at end of file diff --git a/src/axom/mir/docs/sphinx/index.rst b/src/axom/mir/docs/sphinx/index.rst new file mode 100644 index 0000000000..b4b2b257aa --- /dev/null +++ b/src/axom/mir/docs/sphinx/index.rst @@ -0,0 +1,54 @@ +.. ## 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 User Documentation +======================= + +Axom's Material Interface Reconstruction (MIR) component provides algorithms for +reconstructing the interface surfaces between different materials in multimaterial +meshes. The goal of the component is to provide simulation code developers with a +set of algorithms for performing material interface reconstruction with efficient +implementations and an easy-to-use API. + +This component currently provides support for two MIR algorithms: + +.. toctree:: + :maxdepth: 2 + + equi_z + iterative_equi_z + +.. raw:: html + +

Current limitations

+ +.. note:: This MIR component is under active development with many features planned. + +* Improved efficiency for the Equi-Z algorithm implementation is under development. +* Support for GPUs in MIR is planned. +* Support for material interface reconstruction in high-order meshes is planned. + + +.. raw:: html + +

Example Code

+ +In order to access the MIR functionality in this component, import the MIR header file: + +.. literalinclude:: ../../examples/mir_minimal_tutorial.cpp + :start-after: _mir_header_start + :end-before: _mir_header_end + :language: C++ + +The following is all the code needed to generate a mesh with the ``MIRMesh`` class from one of +the default test cases, perform material interface reconstruction using the Equi-Z algorithm +with the ``computeReconstructedInterface()`` function, and then write out the resulting mesh +to a vtk file. + +.. literalinclude:: ../../examples/mir_minimal_tutorial.cpp + :start-after: _mir_main_loop_start + :end-before: _mir_main_loop_end + :language: C++ \ No newline at end of file diff --git a/src/axom/mir/docs/sphinx/iterative_equi_z.rst b/src/axom/mir/docs/sphinx/iterative_equi_z.rst new file mode 100644 index 0000000000..f067cd896b --- /dev/null +++ b/src/axom/mir/docs/sphinx/iterative_equi_z.rst @@ -0,0 +1,24 @@ +.. ## 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) + +========================== +Iterative Equi-Z Algorithm +========================== + +The iterative Equi-Z algorithm provided here is an implementation of the algorithm described +in the 2010 paper "Visualization and Analysis-Oriented Reconstruction of Material Interfaces" +by Meredith and Childs. + +Functionality for this algorithm is provided by ``computeReconstructedInterfaceIterative()``. +It takes in references to two meshes (the mesh to be processed and the mesh that +will contain the results of the material interface reconstruction), an integer denoting the +number of iterations to run for, and a double denoting the percentage difference to modify +the volume fractions by at each iteration. + +This algorithm performs reconstruction using the Equi-Z algorithm, but after finishing recontruction, +calculates a percent difference between the resulting element volume fractions and the desired element volume fractions. +It modifies the original element volume fractions by these values and does another reconstruction with +the Equi-Z algorithm. Since convergence to an optimal value is not guaranteed, it does this for a set number +of iterations that the user specifies. \ No newline at end of file diff --git a/src/axom/mir/docs/sphinx/quad_cases.png b/src/axom/mir/docs/sphinx/quad_cases.png new file mode 100644 index 0000000000..aa6de0cc94 Binary files /dev/null and b/src/axom/mir/docs/sphinx/quad_cases.png differ diff --git a/src/axom/mir/docs/sphinx/triangle_cases.png b/src/axom/mir/docs/sphinx/triangle_cases.png new file mode 100644 index 0000000000..f32e57e941 Binary files /dev/null and b/src/axom/mir/docs/sphinx/triangle_cases.png differ diff --git a/src/axom/mir/docs/sphinx/zoo.png b/src/axom/mir/docs/sphinx/zoo.png new file mode 100644 index 0000000000..49734e028d Binary files /dev/null and b/src/axom/mir/docs/sphinx/zoo.png differ diff --git a/src/axom/mir/examples/CMakeLists.txt b/src/axom/mir/examples/CMakeLists.txt new file mode 100644 index 0000000000..0ec24732d9 --- /dev/null +++ b/src/axom/mir/examples/CMakeLists.txt @@ -0,0 +1,30 @@ +# 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) + +set( mir_examples + mir_tutorial_simple.cpp + mir_minimal_tutorial.cpp + mir_concentric_circles.cpp + ) + +set( mir_example_dependencies + core + slic + mir + ) + +foreach( example ${mir_examples} ) + + get_filename_component( example_name ${example} NAME_WE ) + + blt_add_executable( + NAME ${example_name}_ex + SOURCES ${example} + OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} + DEPENDS_ON ${mir_example_dependencies} + FOLDER axom/mir/examples + ) + +endforeach() diff --git a/src/axom/mir/examples/mir_concentric_circles.cpp b/src/axom/mir/examples/mir_concentric_circles.cpp new file mode 100644 index 0000000000..1455d16f81 --- /dev/null +++ b/src/axom/mir/examples/mir_concentric_circles.cpp @@ -0,0 +1,71 @@ +// 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 "axom/core.hpp" // for axom macros +#include "axom/slic.hpp" +#include "axom/mir.hpp" // for Mir classes & functions + +#include + +// namespace aliases +namespace mir = axom::mir; + +//-------------------------------------------------------------------------------- + +std::string usageString() +{ + return "Args are "; +} + +int main( int argc, char** argv ) +{ + axom::slic::UnitTestLogger logger; // create & initialize test logger + axom::slic::setLoggingMsgLevel( axom::slic::message::Info ); + + if (argc != 4) + { + SLIC_WARNING("Incorrect number of args. " << usageString() ); + return 1; + } + + try + { + // Parse the command line arguments + int gridSize = std::stoi(argv[1]); + int numCircles = std::stoi(argv[2]); + std::string outputFilePath = std::string(argv[3]); + + // Initialize a mesh for testing MIR + auto timer = axom::utilities::Timer(true); + mir::MeshTester tester; + mir::MIRMesh testMesh = tester.initTestCaseFive(gridSize, numCircles); + timer.stop(); + SLIC_INFO("Mesh init time: " << timer.elapsedTimeInMilliSec() << " ms."); + + // Begin material interface reconstruction + timer.start(); + mir::InterfaceReconstructor reconstructor; + mir::MIRMesh processedMesh; + reconstructor.computeReconstructedInterface(testMesh, processedMesh); + timer.stop(); + SLIC_INFO("Material interface reconstruction time: " + << timer.elapsedTimeInMilliSec() << " ms."); + + // Output results + processedMesh.writeMeshToFile(outputFilePath + "outputConcentricCircles.vtk"); + + return 0; + } + catch (std::invalid_argument const &e) + { + SLIC_WARNING("Bad input. " << usageString() ); + return 1; + } + catch (std::out_of_range const &e) + { + SLIC_WARNING("Integer overflow. " << usageString() ); + return 1; + } +} diff --git a/src/axom/mir/examples/mir_minimal_tutorial.cpp b/src/axom/mir/examples/mir_minimal_tutorial.cpp new file mode 100644 index 0000000000..c618928475 --- /dev/null +++ b/src/axom/mir/examples/mir_minimal_tutorial.cpp @@ -0,0 +1,41 @@ +// 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 "axom/core.hpp" +#include "axom/slic.hpp" +#include "axom/slam.hpp" + +// _mir_header_start +#include "axom/mir.hpp" +// _mir_header_end + +#include + +// namespace aliases +namespace mir = axom::mir; +namespace fs = axom::utilities::filesystem; + +/*! + * \brief Tutorial main showing how to initialize test cases and perform mir. + */ +int main(int argc, char** argv) +{ + // _mir_main_loop_start + // Initialize the mesh + mir::MIRMesh testMesh; + mir::MeshTester tester; + testMesh = tester.initTestCaseOne(); + + // Perform material interface reconstruction + mir::MIRMesh processedMesh; + mir::InterfaceReconstructor reconstructor; + reconstructor.computeReconstructedInterface( testMesh, processedMesh ); + + // Output the results to a .vtk file + std::string outputDirectory = fs::joinPath(AXOM_BIN_DIR, "mir_examples/processedMesh.vtk"); + processedMesh.writeMeshToFile( outputDirectory ); + // _mir_main_loop_end + return 0; +} \ No newline at end of file diff --git a/src/axom/mir/examples/mir_tutorial_simple.cpp b/src/axom/mir/examples/mir_tutorial_simple.cpp new file mode 100644 index 0000000000..2de48c4be2 --- /dev/null +++ b/src/axom/mir/examples/mir_tutorial_simple.cpp @@ -0,0 +1,250 @@ +// 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 "axom/core.hpp" +#include "axom/slic.hpp" +#include "axom/slam.hpp" +#include "axom/mir.hpp" + +#include + +// namespace aliases +namespace numerics = axom::numerics; +namespace slam = axom::slam; +namespace mir = axom::mir; +namespace fs = axom::utilities::filesystem; + +//-------------------------------------------------------------------------------- + +enum InputStatus +{ + SUCCESS, + INVALID, + SHOWHELP +}; + +struct Input +{ + int m_test_case; // valid values 1,2,3,4,5 + bool m_should_iterate; + int m_iter_count; + double m_iter_percent; + bool m_verbose; + std::string m_output_dir; + InputStatus m_status; + + Input() + : m_test_case(2) + , m_should_iterate(false) + , m_iter_count(100) + , m_iter_percent(.3) + , m_verbose(false) + , m_status(SUCCESS) + { + m_output_dir = fs::joinPath(AXOM_BIN_DIR, "mir_examples"); + } + + Input(int argc, char** argv) : Input() + { + for (int i = 1 ; i < argc ; /* increment i in loop */) + { + std::string arg = argv[i]; + if (arg == "--test-case") + { + m_test_case = std::stoi(argv[++i]); + } + else if (arg == "--output-dir") + { + m_output_dir = argv[++i]; + } + else if (arg == "--iter-count") + { + m_iter_count = std::stoi(argv[++i]); + m_should_iterate = true; + } + else if (arg == "--iter-percent") + { + m_iter_percent = std::stod(argv[++i]); + m_should_iterate = true; + } + else if (arg == "--verbose") + { + m_verbose = true; + } + else // help or unknown parameter + { + if(arg != "--help" && arg != "-h") + { + SLIC_WARNING("Unrecognized parameter: " << arg); + m_status = INVALID; + } + else + { + m_status = SHOWHELP; + } + return; + } + ++i; + } + + checkTestCase(); + checkOutputDir(); + checkIterationParams(); + } + + + bool shouldIterate() const { return m_should_iterate; } + int numIterations() const { return m_iter_count; } + double iterPercentage() const { return m_iter_percent; } + + + void showhelp() + { + std::cout + << "Argument usage:" + "\n --help Show this help message." + "\n --test-case N Mesh test case. Default N = 2." + "\n Valid values {1,2,3,4,5}" + "\n --output-dir dir Directory for output mesh" + "\n Default is: '${AXOM_BIN_DIR}/mir_examples'" + "\n --iter-count N Number of iterations for iterative algorithm" + "\n Defaults to 100." + "\n Setting a value triggers iterative algorithm" + "\n --iter-percent D Volume diff percentage for iterative algorithm" + "\n Must be between 0 and 1. Defaults to 0.3" + "\n Setting a value triggers iterative algorithm" + "\n --verbose Increases verbosity of output" + << std::endl << std::endl; + }; + +private: + void checkTestCase() + { + if(m_test_case < 1 || m_test_case > 5) + { + m_status = INVALID; + SLIC_WARNING("Invalid test case " << m_test_case); + } + } + + void checkOutputDir() + { + if(! fs::pathExists(m_output_dir)) + { + fs::makeDirsForPath(m_output_dir); + } + } + + void checkIterationParams() + { + if(m_should_iterate) + { + if(m_iter_count < 1) + { + m_status = INVALID; + SLIC_WARNING("Invalid iteration count " << m_iter_count); + } + + if(m_iter_percent <= 0. || m_iter_percent > 1.) + { + m_status = INVALID; + SLIC_WARNING("Invalid iteration percentage " << m_iter_percent); + } + } + } + +}; + +/*! + * \brief Tutorial main showing how to initialize test cases and perform mir. + */ +int main(int argc, char** argv) +{ + axom::slic::UnitTestLogger logger; // create & initialize test logger + axom::slic::setLoggingMsgLevel( axom::slic::message::Info ); + + // Parse arguments + Input params(argc, argv); + + if (params.m_status != SUCCESS) + { + if (params.m_status == SHOWHELP) + { + params.showhelp(); + return 0; + } + else if (params.m_status == INVALID) + { + params.showhelp(); + return 1; + } + } + + mir::MIRMesh testMesh; + mir::MeshTester tester; + + auto timer = axom::utilities::Timer(true); + + switch(params.m_test_case) + { + case 1: testMesh = tester.initTestCaseOne(); break; + case 2: testMesh = tester.initTestCaseTwo(); break; + case 3: testMesh = tester.initTestCaseThree(); break; + case 4: testMesh = tester.initTestCaseFour(); break; + case 5: testMesh = tester.initTestCaseFive(25, 12); break; + } + + timer.stop(); + SLIC_INFO( "Mesh init time: " << timer.elapsedTimeInMilliSec() << " ms."); + + SLIC_INFO("Test mesh is " << (testMesh.isValid(true) ? "" : " NOT") << " valid."); + + if(params.m_verbose) + { + SLIC_INFO("Initial mesh:"); + testMesh.print(); + } + + // Begin material interface reconstruction + timer.start(); + + mir::MIRMesh processedMesh; + mir::InterfaceReconstructor reconstructor; + + if(! params.shouldIterate()) // Process once, with original Meredith algorithm + { + reconstructor.computeReconstructedInterface(testMesh, processedMesh); + } + else // use iterative algorithm + { + int n = params.numIterations(); + double p = params.iterPercentage(); + + reconstructor.computeReconstructedInterfaceIterative(testMesh, n, p, processedMesh); + } + + timer.stop(); + SLIC_INFO( "Reconstruction time: " << timer.elapsedTimeInMilliSec() << " ms."); + + // Output results + std::string dir = params.m_output_dir; + processedMesh.writeMeshToFile( fs::joinPath(dir, "processedMesh.vtk") ); + + using VolFracs = std::vector >; + timer.start(); + VolFracs materialVolumeFractionsElement = processedMesh.computeOriginalElementVolumeFractions(); + timer.stop(); + SLIC_INFO( "Computing volumes took: " << timer.elapsedTimeInMilliSec() << " ms."); + + if(params.m_verbose) + { + SLIC_INFO("Final mesh:"); + processedMesh.print(); + } + + return 0; +} + +//-------------------------------------------------------------------------------- diff --git a/src/axom/mir/tests/CMakeLists.txt b/src/axom/mir/tests/CMakeLists.txt new file mode 100644 index 0000000000..d4b750967b --- /dev/null +++ b/src/axom/mir/tests/CMakeLists.txt @@ -0,0 +1,42 @@ +# 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 unit tests +#------------------------------------------------------------------------------ + + +#------------------------------------------------------------------------------ +# Specify list of tests +#------------------------------------------------------------------------------ + +set(gtest_mir_tests + mir_smoke.cpp + mir_mesh.cpp + mir_interface_reconstructor.cpp + mir_cell_clipper.cpp + mir_utilities.cpp + mir_cell_generator.cpp + ) + +set(mir_tests_depends_on + slic + mir + gtest) + +#------------------------------------------------------------------------------ +# Add gtest based tests +#------------------------------------------------------------------------------ +foreach(test ${gtest_mir_tests}) + get_filename_component( test_name ${test} NAME_WE ) + blt_add_executable( NAME ${test_name}_test + SOURCES ${test} + OUTPUT_DIR ${TEST_OUTPUT_DIRECTORY} + DEPENDS_ON ${mir_tests_depends_on} + FOLDER axom/mir/tests ) + + blt_add_test( NAME ${test_name} + COMMAND ${test_name}_test ) +endforeach() + diff --git a/src/axom/mir/tests/mir_cell_clipper.cpp b/src/axom/mir/tests/mir_cell_clipper.cpp new file mode 100644 index 0000000000..68f538e3f4 --- /dev/null +++ b/src/axom/mir/tests/mir_cell_clipper.cpp @@ -0,0 +1,451 @@ +// 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) + +#ifndef MIR_CELL_CLIPPER_TEST_H_ +#define MIR_CELL_CLIPPER_TEST_H_ + +#include "gtest/gtest.h" + +#include "axom/slic.hpp" +#include "axom/mir.hpp" + +using namespace axom; + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, triangle_case_zero) +{ + mir::Shape shape = mir::Shape::Triangle; + std::vector matOneVF = {0.0, 0.0, 0.0}; + std::vector matTwoVF = {1.0, 1.0, 1.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + + EXPECT_EQ( 0 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, triangle_case_three) +{ + mir::Shape shape = mir::Shape::Triangle; + std::vector matOneVF = {0.0, 1.0, 1.0}; + std::vector matTwoVF = {1.0, 0.0, 0.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + EXPECT_EQ( 3 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, triangle_case_four) +{ + mir::Shape shape = mir::Shape::Triangle; + std::vector matOneVF = {1.0, 0.0, 0.0}; + std::vector matTwoVF = {0.0, 1.0, 1.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + + EXPECT_EQ( 4 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, triangle_case_seven) +{ + mir::Shape shape = mir::Shape::Triangle; + std::vector matOneVF = {1.0, 1.0, 1.0}; + std::vector matTwoVF = {0.0, 0.0, 0.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + + EXPECT_EQ( 7 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, quad_case_zero) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {0.0, 0.0, 0.0, 0.0}; + std::vector matTwoVF = {1.0, 1.0, 1.0, 1.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + EXPECT_EQ( 0 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, quad_case_one) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {0.0, 0.0, 0.0, 1.0}; + std::vector matTwoVF = {1.0, 1.0, 1.0, 0.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + EXPECT_EQ( 1 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, quad_case_two) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {0.0, 0.0, 1.0, 0.0}; + std::vector matTwoVF = {1.0, 1.0, 0.0, 1.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + EXPECT_EQ( 2 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, quad_case_three) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {0.0, 0.0, 1.0, 1.0}; + std::vector matTwoVF = {1.0, 1.0, 0.0, 0.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + EXPECT_EQ( 3 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, quad_case_five) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {0.0, 1.0, 0.0, 1.0}; + std::vector matTwoVF = {1.0, 0.0, 1.0, 0.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + EXPECT_EQ( 5 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, quad_case_ten) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {1.0, 0.0, 1.0, 0.0}; + std::vector matTwoVF = {0.0, 1.0, 0.0, 1.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + EXPECT_EQ( 10 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_case, quad_case_fifteen) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {1.0, 1.0, 1.0, 1.0}; + std::vector matTwoVF = {0.0, 0.0, 0.0, 0.0}; + + mir::CellClipper clipper; + unsigned int clippingCase = clipper.determineClippingCase(shape, matOneVF, matTwoVF); + EXPECT_EQ( 15 , clippingCase); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_place, clip_edge_when_mat_one_dominates) +{ + axom::float64 vfMatOneVertexOne = 0.0; + axom::float64 vfMatTwoVertexOne = 1.0; + + axom::float64 vfMatOneVertexTwo = 0.0; + axom::float64 vfMatTwoVertexTwo = 1.0; + + mir::CellClipper clipper; + axom::float64 tValue = clipper.computeTValueOnEdge(vfMatOneVertexOne, vfMatTwoVertexOne, vfMatOneVertexTwo, vfMatTwoVertexTwo); + + EXPECT_DOUBLE_EQ( tValue, 0.0 ); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_place, clip_edge_when_mat_two_dominates) +{ + axom::float64 vfMatOneVertexOne = 1.0; + axom::float64 vfMatTwoVertexOne = 0.0; + + axom::float64 vfMatOneVertexTwo = 1.0; + axom::float64 vfMatTwoVertexTwo = 0.0; + + mir::CellClipper clipper; + axom::float64 tValue = clipper.computeTValueOnEdge(vfMatOneVertexOne, vfMatTwoVertexOne, vfMatOneVertexTwo, vfMatTwoVertexTwo); + + EXPECT_DOUBLE_EQ( tValue, 0.0 ); +} + +// //---------------------------------------------------------------------- + + +TEST(mir_clipping_place, clip_edge_in_middle) +{ + axom::float64 vfMatOneVertexOne = 1.0; + axom::float64 vfMatTwoVertexOne = 0.0; + + axom::float64 vfMatOneVertexTwo = 0.0; + axom::float64 vfMatTwoVertexTwo = 1.0; + + mir::CellClipper clipper; + axom::float64 tValue = clipper.computeTValueOnEdge(vfMatOneVertexOne, vfMatTwoVertexOne, vfMatOneVertexTwo, vfMatTwoVertexTwo); + + EXPECT_DOUBLE_EQ( tValue, 0.5 ); +} + +// //---------------------------------------------------------------------- + +TEST(mir_clipping_place, clip_edge_with_one_null_material) +{ + axom::float64 vfMatOneVertexOne = -1.0; + axom::float64 vfMatTwoVertexOne = 0.0; + + axom::float64 vfMatOneVertexTwo = -1.0; + axom::float64 vfMatTwoVertexTwo = 0.0; + + mir::CellClipper clipper; + axom::float64 tValue = clipper.computeTValueOnEdge(vfMatOneVertexOne, vfMatTwoVertexOne, vfMatOneVertexTwo, vfMatTwoVertexTwo); + + EXPECT_DOUBLE_EQ( tValue, 0.0 ); +} + +// //---------------------------------------------------------------------- + +TEST(mir_clipping_place, clip_edge_with_two_null_material) +{ + axom::float64 vfMatOneVertexOne = -1.0; + axom::float64 vfMatTwoVertexOne = -1.0; + + axom::float64 vfMatOneVertexTwo = -1.0; + axom::float64 vfMatTwoVertexTwo = -1.0; + + mir::CellClipper clipper; + axom::float64 tValue = clipper.computeTValueOnEdge(vfMatOneVertexOne, vfMatTwoVertexOne, vfMatOneVertexTwo, vfMatTwoVertexTwo); + + EXPECT_DOUBLE_EQ( tValue, 0.0 ); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_cell_and_vertex_output, clip_quad_case_zero) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {0.0, 0.0, 0.0, 0.0}; + std::vector matTwoVF = {1.0, 1.0, 1.0, 1.0}; + + std::vector > vertexVF; + vertexVF.push_back(matOneVF); + vertexVF.push_back(matTwoVF); + + std::map > newElements; + std::map > newVertices; + axom::float64 tValues[8] = { 0 }; + + mir::CellClipper clipper; + clipper.computeClippingPoints(shape, vertexVF, newElements, newVertices, tValues); + + EXPECT_EQ( 1, newElements.size() ); + EXPECT_EQ( 4, newVertices.size() ); + + EXPECT_EQ( 0, newElements[0][0] ); + EXPECT_EQ( 1, newElements[0][1] ); + EXPECT_EQ( 2, newElements[0][2] ); + EXPECT_EQ( 3, newElements[0][3] ); + + EXPECT_EQ( 0, newVertices[0][0] ); + EXPECT_EQ( 0, newVertices[1][0] ); + EXPECT_EQ( 0, newVertices[2][0] ); + EXPECT_EQ( 0, newVertices[3][0] ); + + for (int i = 0; i < 8; ++i) + { + EXPECT_DOUBLE_EQ( 0, tValues[i] ); + } +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_cell_and_vertex_output, clip_quad_case_one) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {0.0, 0.0, 0.0, 1.0}; + std::vector matTwoVF = {1.0, 1.0, 1.0, 0.0}; + + std::vector > vertexVF; + vertexVF.push_back(matOneVF); + vertexVF.push_back(matTwoVF); + + std::map > newElements; + std::map > newVertices; + axom::float64 tValues[8] = { 0 }; + + mir::CellClipper clipper; + clipper.computeClippingPoints(shape, vertexVF, newElements, newVertices, tValues); + + EXPECT_EQ( 3, newElements.size() ); + EXPECT_EQ( 6, newVertices.size() ); + + // Check the first element + EXPECT_EQ( 3, newElements[0][0] ); + EXPECT_EQ( 7, newElements[0][1] ); + EXPECT_EQ( 6, newElements[0][2] ); + + // Check the second element + EXPECT_EQ( 7, newElements[1][0] ); + EXPECT_EQ( 0, newElements[1][1] ); + EXPECT_EQ( 2, newElements[1][2] ); + EXPECT_EQ( 6, newElements[1][3] ); + + // Check the third element + EXPECT_EQ( 0, newElements[2][0] ); + EXPECT_EQ( 1, newElements[2][1] ); + EXPECT_EQ( 2, newElements[2][2] ); + + // Check each vertex's associated elements + EXPECT_EQ( 1, newVertices[0][0] ); + EXPECT_EQ( 2, newVertices[0][1] ); + + EXPECT_EQ( 2, newVertices[1][0] ); + + EXPECT_EQ( 1, newVertices[2][0] ); + EXPECT_EQ( 2, newVertices[2][1] ); + + EXPECT_EQ( 0, newVertices[3][0] ); + + EXPECT_EQ( 0, newVertices[6][0] ); + EXPECT_EQ( 1, newVertices[6][1] ); + + EXPECT_EQ( 0, newVertices[7][0] ); + EXPECT_EQ( 1, newVertices[7][1] ); +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_cell_and_vertex_output, clip_quad_case_three) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {0.0, 0.0, 1.0, 1.0}; + std::vector matTwoVF = {1.0, 1.0, 0.0, 0.0}; + + std::vector > vertexVF; + vertexVF.push_back(matOneVF); + vertexVF.push_back(matTwoVF); + + std::map > newElements; + std::map > newVertices; + axom::float64 tValues[8] = { 0 }; + + mir::CellClipper clipper; + clipper.computeClippingPoints(shape, vertexVF, newElements, newVertices, tValues); + + EXPECT_EQ( 2, newElements.size() ); + EXPECT_EQ( 6, newVertices.size() ); + + EXPECT_EQ( 0, newElements[0][0] ); + EXPECT_EQ( 1, newElements[0][1] ); + EXPECT_EQ( 5, newElements[0][2] ); + EXPECT_EQ( 7, newElements[0][3] ); + + EXPECT_EQ( 7, newElements[1][0] ); + EXPECT_EQ( 5, newElements[1][1] ); + EXPECT_EQ( 2, newElements[1][2] ); + EXPECT_EQ( 3, newElements[1][3] ); + + EXPECT_EQ( 0, newVertices[0][0] ); + EXPECT_EQ( 0, newVertices[1][0] ); + EXPECT_EQ( 1, newVertices[2][0] ); + EXPECT_EQ( 1, newVertices[3][0] ); + EXPECT_EQ( 0, newVertices[5][0] ); + EXPECT_EQ( 1, newVertices[5][1] ); + EXPECT_EQ( 0, newVertices[7][0] ); + EXPECT_EQ( 1, newVertices[7][1] ); + +} + +//---------------------------------------------------------------------- + +TEST(mir_clipping_cell_and_vertex_output, clip_quad_case_five) +{ + mir::Shape shape = mir::Shape::Quad; + std::vector matOneVF = {0.0, 1.0, 0.0, 1.0}; + std::vector matTwoVF = {1.0, 0.0, 1.0, 0.0}; + + std::vector > vertexVF; + vertexVF.push_back(matOneVF); + vertexVF.push_back(matTwoVF); + + std::map > newElements; + std::map > newVertices; + axom::float64 tValues[8] = { 0 }; + + mir::CellClipper clipper; + clipper.computeClippingPoints(shape, vertexVF, newElements, newVertices, tValues); + + EXPECT_EQ( 4, newElements.size() ); + EXPECT_EQ( 8, newVertices.size() ); + + EXPECT_EQ( 4, newElements[0][0] ); + EXPECT_EQ( 1, newElements[0][1] ); + EXPECT_EQ( 5, newElements[0][2] ); + + EXPECT_EQ( 0, newElements[1][0] ); + EXPECT_EQ( 4, newElements[1][1] ); + EXPECT_EQ( 5, newElements[1][2] ); + EXPECT_EQ( 2, newElements[1][3] ); + + EXPECT_EQ( 0, newElements[2][0] ); + EXPECT_EQ( 2, newElements[2][1] ); + EXPECT_EQ( 6, newElements[2][2] ); + EXPECT_EQ( 7, newElements[2][3] ); + + EXPECT_EQ( 7, newElements[3][0] ); + EXPECT_EQ( 6, newElements[3][1] ); + EXPECT_EQ( 3, newElements[3][2] ); + + EXPECT_EQ( 1, newVertices[0][0] ); + EXPECT_EQ( 2, newVertices[0][1] ); + EXPECT_EQ( 0, newVertices[1][0] ); + EXPECT_EQ( 1, newVertices[2][0] ); + EXPECT_EQ( 2, newVertices[2][1] ); + EXPECT_EQ( 3, newVertices[3][0] ); + EXPECT_EQ( 0, newVertices[4][0] ); + EXPECT_EQ( 1, newVertices[4][1] ); + EXPECT_EQ( 0, newVertices[5][0] ); + EXPECT_EQ( 1, newVertices[5][1] ); + EXPECT_EQ( 2, newVertices[6][0] ); + EXPECT_EQ( 3, newVertices[6][1] ); + EXPECT_EQ( 2, newVertices[7][0] ); + EXPECT_EQ( 3, newVertices[7][1] ); + +} + +//---------------------------------------------------------------------- + +int main(int argc, char* argv[]) +{ + int result = 0; + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::UnitTestLogger logger; // create & initialize test logger, + + result = RUN_ALL_TESTS(); + return result; +} + + +#endif // MIR_CELL_CLIPPER_TEST_H_ diff --git a/src/axom/mir/tests/mir_cell_generator.cpp b/src/axom/mir/tests/mir_cell_generator.cpp new file mode 100644 index 0000000000..9505eae633 --- /dev/null +++ b/src/axom/mir/tests/mir_cell_generator.cpp @@ -0,0 +1,231 @@ +// 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) + +#ifndef MIR_CELL_GENERATOR_TEST_H_ +#define MIR_CELL_GENERATOR_TEST_H_ + +#include "gtest/gtest.h" + +#include "axom/slic.hpp" +#include "axom/mir.hpp" + +using namespace axom; + +//---------------------------------------------------------------------- + +TEST(mir_cell_generator, generate_quad_topology) +{ + EXPECT_EQ(true, 1); + // this function generates the evInds, evBegins, ... etc given the map of elements -> verts, and verts -> elements + std::map > elementMap; + elementMap[0] = { 4, 1, 5 }; + elementMap[1] = { 0, 4, 5, 2 }; + elementMap[2] = { 0, 2, 6, 7 }; + elementMap[3] = { 7, 6, 3 }; + + std::map > vertexMap; + vertexMap[0] = { 1, 2 }; + vertexMap[1] = { 0 }; + vertexMap[2] = { 1, 2 }; + vertexMap[3] = { 3 }; + vertexMap[4] = { 0, 1 }; + vertexMap[5] = { 0, 1 }; + vertexMap[6] = { 2, 3 }; + vertexMap[7] = { 2, 3 }; + + mir::CellData cellData; + mir::CellGenerator generator; + generator.generateTopologyData(elementMap, vertexMap, cellData); + + EXPECT_EQ( 4, cellData.m_topology.m_evInds[0] ); + EXPECT_EQ( 1, cellData.m_topology.m_evInds[1] ); + EXPECT_EQ( 5, cellData.m_topology.m_evInds[2] ); + EXPECT_EQ( 0, cellData.m_topology.m_evInds[3] ); + EXPECT_EQ( 4, cellData.m_topology.m_evInds[4] ); + EXPECT_EQ( 5, cellData.m_topology.m_evInds[5] ); + EXPECT_EQ( 2, cellData.m_topology.m_evInds[6] ); + EXPECT_EQ( 0, cellData.m_topology.m_evInds[7] ); + EXPECT_EQ( 2, cellData.m_topology.m_evInds[8] ); + EXPECT_EQ( 6, cellData.m_topology.m_evInds[9] ); + EXPECT_EQ( 7, cellData.m_topology.m_evInds[10] ); + EXPECT_EQ( 7, cellData.m_topology.m_evInds[11] ); + EXPECT_EQ( 6, cellData.m_topology.m_evInds[12] ); + EXPECT_EQ( 3, cellData.m_topology.m_evInds[13] ); + + EXPECT_EQ( 1, cellData.m_topology.m_veInds[0] ); + EXPECT_EQ( 2, cellData.m_topology.m_veInds[1] ); + EXPECT_EQ( 0, cellData.m_topology.m_veInds[2] ); + EXPECT_EQ( 1, cellData.m_topology.m_veInds[3] ); + EXPECT_EQ( 2, cellData.m_topology.m_veInds[4] ); + EXPECT_EQ( 3, cellData.m_topology.m_veInds[5] ); + EXPECT_EQ( 0, cellData.m_topology.m_veInds[6] ); + EXPECT_EQ( 1, cellData.m_topology.m_veInds[7] ); + EXPECT_EQ( 0, cellData.m_topology.m_veInds[8] ); + EXPECT_EQ( 1, cellData.m_topology.m_veInds[9] ); + EXPECT_EQ( 2, cellData.m_topology.m_veInds[10] ); + EXPECT_EQ( 3, cellData.m_topology.m_veInds[11] ); + EXPECT_EQ( 2, cellData.m_topology.m_veInds[12] ); + EXPECT_EQ( 3, cellData.m_topology.m_veInds[13] ); + + EXPECT_EQ( 0, cellData.m_topology.m_evBegins[0] ); + EXPECT_EQ( 3, cellData.m_topology.m_evBegins[1] ); + EXPECT_EQ( 7, cellData.m_topology.m_evBegins[2] ); + EXPECT_EQ( 11, cellData.m_topology.m_evBegins[3] ); + + EXPECT_EQ( 0, cellData.m_topology.m_veBegins[0] ); + EXPECT_EQ( 2, cellData.m_topology.m_veBegins[1] ); + EXPECT_EQ( 3, cellData.m_topology.m_veBegins[2] ); + EXPECT_EQ( 5, cellData.m_topology.m_veBegins[3] ); + EXPECT_EQ( 6, cellData.m_topology.m_veBegins[4] ); + EXPECT_EQ( 8, cellData.m_topology.m_veBegins[5] ); + EXPECT_EQ( 10, cellData.m_topology.m_veBegins[6] ); + EXPECT_EQ( 12, cellData.m_topology.m_veBegins[7] ); +} + +//---------------------------------------------------------------------- + +TEST(mir_cell_generator, generate_vertex_positions) +{ + mir::Shape shapeType = mir::Shape::Quad; + + std::map > vertexMap; + vertexMap[0] = { 1, 2 }; + vertexMap[1] = { 0 }; + vertexMap[2] = { 1, 2 }; + vertexMap[3] = { 3 }; + vertexMap[4] = { 0, 1 }; + vertexMap[5] = { 0, 1 }; + vertexMap[6] = { 2, 3 }; + vertexMap[7] = { 2, 3 }; + + std::vector originalVertexPositions; + originalVertexPositions.push_back( mir::Point2::make_point(0.0, 1.0) ); + originalVertexPositions.push_back( mir::Point2::make_point(0.0, 0.0) ); + originalVertexPositions.push_back( mir::Point2::make_point(1.0, 0.0) ); + originalVertexPositions.push_back( mir::Point2::make_point(1.0, 1.0) ); + + axom::float64 tValues[8] = { 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5 }; + + mir::CellData cellData; + mir::CellGenerator cellGenerator; + cellGenerator.generateVertexPositions(shapeType, vertexMap, originalVertexPositions, tValues, cellData); + + const auto& positions = cellData.m_mapData.m_vertexPositions; + EXPECT_NEAR( positions[0][0], 0.0, 0.00001 ); + EXPECT_NEAR( positions[0][1], 1.0, 0.00001 ); + + EXPECT_NEAR( positions[1][0], 0.0, 0.00001 ); + EXPECT_NEAR( positions[1][1], 0.0, 0.00001 ); + + EXPECT_NEAR( positions[2][0], 1.0, 0.00001 ); + EXPECT_NEAR( positions[2][1], 0.0, 0.00001 ); + + EXPECT_NEAR( positions[3][0], 1.0, 0.00001 ); + EXPECT_NEAR( positions[3][1], 1.0, 0.00001 ); + + EXPECT_NEAR( positions[4][0], 0.0, 0.00001 ); + EXPECT_NEAR( positions[4][1], 0.5, 0.00001 ); + + EXPECT_NEAR( positions[5][0], 0.5, 0.00001 ); + EXPECT_NEAR( positions[5][1], 0.0, 0.00001 ); + + EXPECT_NEAR( positions[6][0], 1.0, 0.00001 ); + EXPECT_NEAR( positions[6][1], 0.5, 0.00001 ); + + EXPECT_NEAR( positions[7][0], 0.5, 0.00001 ); + EXPECT_NEAR( positions[7][1], 1.0, 0.00001 ); +} + +//---------------------------------------------------------------------- + +TEST(mir_cell_generator, generate_vertex_volume_fractions) +{ + mir::Shape shapeType = mir::Shape::Quad; + + std::map > vertexMap; + vertexMap[0] = { 1, 2 }; + vertexMap[1] = { 0 }; + vertexMap[2] = { 1, 2 }; + vertexMap[3] = { 3 }; + vertexMap[4] = { 0, 1 }; + vertexMap[5] = { 0, 1 }; + vertexMap[6] = { 2, 3 }; + vertexMap[7] = { 2, 3 }; + + std::vector > originalVertexVF(2); + originalVertexVF[0] = { 0.0, 0.33, 0.67, 1.0 }; + originalVertexVF[1] = { 1.0, 0.67, 0.33, 0.0 }; + + axom::float64 tValues[8] = { 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5 }; + + mir::CellData cellData; + mir::CellGenerator cellGenerator; + cellGenerator.generateVertexVolumeFractions(shapeType, vertexMap, originalVertexVF, tValues, cellData); + + + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[0][0], 0.0, 0.00001 ); + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[1][0], 1.0, 0.00001 ); + + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[0][1], 0.33, 0.00001 ); + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[1][1], 0.67, 0.00001 ); + + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[0][2], 0.67, 0.00001 ); + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[1][2], 0.33, 0.00001 ); + + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[0][3], 1.0, 0.00001 ); + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[1][3], 0.0, 0.00001 ); + + + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[0][4], 0.165, 0.00001 ); + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[1][4], 0.835, 0.00001 ); + + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[0][5], 0.5, 0.00001 ); + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[1][5], 0.5, 0.00001 ); + + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[0][6], 0.835, 0.00001 ); + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[1][6], 0.165, 0.00001 ); + + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[0][7], 0.5, 0.00001 ); + EXPECT_NEAR( cellData.m_mapData.m_vertexVolumeFractions[1][7], 0.5, 0.00001 ); +} + +//---------------------------------------------------------------------- + +TEST(mir_cell_generator, determine_clean_cell_material) +{ + mir::Shape shapeType = mir::Shape::Quad; + + std::vector vertexIDs = { 0, 1, 5, 7 }; + + int matOne = 0; + int matTwo = 1; + + std::vector > originalVertexVF(2); + originalVertexVF[0] = { 0.0, 0.33, 0.67, 1.0 }; + originalVertexVF[1] = { 1.0, 0.67, 0.33, 0.0 }; + + mir::CellData cellData; + mir::CellGenerator cellGenerator; + + int dominantMaterial = cellGenerator.determineCleanCellMaterial(shapeType, vertexIDs, matOne, matTwo, originalVertexVF); + + EXPECT_EQ( dominantMaterial, 1 ); +} + +//---------------------------------------------------------------------- + +int main(int argc, char* argv[]) +{ + int result = 0; + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::UnitTestLogger logger; // create & initialize test logger, + + result = RUN_ALL_TESTS(); + return result; +} + + +#endif // MIR_CELL_GENERATOR_TEST_H_ diff --git a/src/axom/mir/tests/mir_interface_reconstructor.cpp b/src/axom/mir/tests/mir_interface_reconstructor.cpp new file mode 100644 index 0000000000..875c7365df --- /dev/null +++ b/src/axom/mir/tests/mir_interface_reconstructor.cpp @@ -0,0 +1,30 @@ +// 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) + +#ifndef MIR_INTERFACE_RECONSTRUCTOR_TEST_H_ +#define MIR_INTERFACE_RECONSTRUCTOR_TEST_H_ + +#include "gtest/gtest.h" + +#include "axom/slic.hpp" +#include "axom/mir.hpp" + +using namespace axom; + +//---------------------------------------------------------------------- + +int main(int argc, char* argv[]) +{ + int result = 0; + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::UnitTestLogger logger; // create & initialize test logger, + + result = RUN_ALL_TESTS(); + return result; +} + + +#endif // MIR_INTERFACE_RECONSTRUCTOR_TEST_H_ diff --git a/src/axom/mir/tests/mir_mesh.cpp b/src/axom/mir/tests/mir_mesh.cpp new file mode 100644 index 0000000000..0d9797e16d --- /dev/null +++ b/src/axom/mir/tests/mir_mesh.cpp @@ -0,0 +1,229 @@ +// 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) + +#ifndef MIR_MESH_TEST_H_ +#define MIR_MESH_TEST_H_ + +#include "gtest/gtest.h" + +#include "axom/slic.hpp" +#include "axom/mir.hpp" + + +namespace mir = axom::mir; + +class MirMeshTest : public ::testing::Test +{ +public: + + template + using Vec = std::vector; + + using IndexVec = Vec; + using VolFracVec = Vec; + using VolumeFractions = Vec; + + enum { GREEN = 0, BLUE = 1 }; + +public: + + mir::MIRMesh& getMesh() { return m_mesh; } + const mir::MIRMesh& getMesh() const { return m_mesh; } + +protected: + void SetUp() override + { + m_verts= mir::VertSet(16); + m_elems= mir::ElemSet(9); + m_nMats = 2; + + // Create the mesh connectivity information + setupTopoData(); + setupMapData(); + setupVolumeFractions(); + + m_mesh.initializeMesh(this->m_verts, this->m_elems, this->m_nMats, + this->m_topoData, this->m_mapData, this->m_volFracs); + } + + // Set up the mesh's topological connectivity + void setupTopoData() + { + m_topoData.m_evInds = { + 0,4,5,1, // elem 0, card 4, start 0 + 1,5,6,2, // elem 1, card 4, start 4 + 2,6,7,3, // elem 2, card 4, start 8 + 4,8,9,5, // elem 3, card 4, start 12 + 5,9,10,6, // elem 4, card 4, start 16 + 6,10,11,7, // elem 5, card 4, start 20 + 8,12,13,9, // elem 6, card 4, start 24 + 9,13,14,10, // elem 7, card 4, start 28 + 10,14,15,11 // elem 8, card 4, start 32, end 36 + }; + + m_topoData.m_evBegins = { + 0,4,8,12,16,20,24,28,32,36 + }; + + m_topoData.m_veInds = { + 0, // vert 0, card 1, start 0 + 0,1, // vert 1, card 2, start 1 + 1,2, // vert 2, card 2, start 3 + 2, // vert 3, card 1, start 5 + 0,3, // vert 4, card 2, start 6 + 0,1,3,4, // vert 5, card 4, start 8 + 1,2,4,5, // vert 6, card 4, start 12 + 2,5, // vert 7, card 2, start 16 + 3,6, // vert 8, card 2, start 18 + 3,4,6,7, // vert 9, card 4, start 20 + 4,5,7,8, // vert 10, card 4, start 24 + 5,8, // vert 11, card 2, start 28 + 6, // vert 12, card 1, start 30 + 6,7, // vert 13, card 2, start 31 + 7,8, // vert 14, card 2, start 33 + 8, // vert 15, card 1, start 35, end 36 + }; + + m_topoData.m_veBegins = { + 0,1,3,5,6,8,12,16,18,20,24,28,30,31,33,35,36 + }; + } + + // Set up the mesh's map data + void setupMapData() + { + m_mapData.m_vertexPositions = { + mir::Point2( 0.0, 3.0 ), + mir::Point2( 1.0, 3.0 ), + mir::Point2( 2.0, 3.0 ), + mir::Point2( 3.0, 3.0 ), + + mir::Point2( 0.0, 2.0 ), + mir::Point2( 1.0, 2.0 ), + mir::Point2( 2.0, 2.0 ), + mir::Point2( 3.0, 2.0 ), + + mir::Point2( 0.0, 1.0 ), + mir::Point2( 1.0, 1.0 ), + mir::Point2( 2.0, 1.0 ), + mir::Point2( 3.0, 1.0 ), + + mir::Point2( 0.0, 0.0 ), + mir::Point2( 1.0, 0.0 ), + mir::Point2( 2.0, 0.0 ), + mir::Point2( 3.0, 0.0 ) + }; + + m_mapData.m_elementDominantMaterials = Vec(m_elems.size(), mir::NULL_MAT); + m_mapData.m_elementParents = { 0, 1, 2, 3, 4, 5, 6, 7, 8 }; + m_mapData.m_shapeTypes = Vec(m_elems.size(), mir::Shape::Quad); + } + + // Set up the mesh's volume fraction data + void setupVolumeFractions() + { + m_volFracs.resize(m_nMats); + m_volFracs[GREEN] = {1.0, 1.0, 1.0, 1.0, 0.5, 0.2, 0.2, 0.0, 0.0}; + m_volFracs[BLUE] = {0.0, 0.0, 0.0, 0.0, 0.5, 0.8, 0.8, 1.0, 1.0}; + } + +protected: + mir::VertSet m_verts; + mir::ElemSet m_elems; + int m_nMats; + + mir::CellTopologyData m_topoData; + mir::CellMapData m_mapData; + mir::CellData m_cellData; + VolumeFractions m_volFracs; + + mir::MIRMesh m_mesh; +}; + + +TEST_F(MirMeshTest,default_ctor) +{ + mir::MIRMesh mesh; + EXPECT_TRUE(mesh.isValid(true)); +} + + +TEST_F(MirMeshTest,initialize) +{ + mir::MIRMesh mesh; + mesh.initializeMesh( + this->m_verts, + this->m_elems, + this->m_nMats, + this->m_topoData, + this->m_mapData, + this->m_volFracs); + + EXPECT_TRUE(mesh.isValid(true)); + + mesh.print(); +} + +TEST_F(MirMeshTest,copy_ctor) +{ + // test copy constructor + { + mir::MIRMesh& mesh = this->getMesh(); + EXPECT_TRUE(mesh.isValid(true)); + + mir::MIRMesh mirCopy(mesh); + EXPECT_TRUE( mirCopy.isValid(true) ); + } + + // test const copy constructor + { + const mir::MIRMesh& cmesh = this->getMesh(); + EXPECT_TRUE(cmesh.isValid(true)); + + // test copy constructor + const mir::MIRMesh mirCopy(cmesh); + EXPECT_TRUE( mirCopy.isValid(true) ); + } +} + +TEST_F(MirMeshTest,copy_assign) +{ + // test copy assignment from a non-const MIRMesh + { + mir::MIRMesh& mesh = this->getMesh(); + EXPECT_TRUE(mesh.isValid(true)); + + mir::MIRMesh mirCopy; + mirCopy = mesh; + EXPECT_TRUE( mirCopy.isValid(true) ); + } + + // test copy assignment from a const MIRMesh + { + const mir::MIRMesh& cmesh = this->getMesh(); + EXPECT_TRUE(cmesh.isValid(true)); + + mir::MIRMesh mirCopy; + mirCopy = cmesh; + EXPECT_TRUE( mirCopy.isValid(true) ); + } +} + + +//------------------------------------------------------------------------------ + +int main(int argc, char* argv[]) +{ + int result = 0; + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::UnitTestLogger logger; // create & initialize test logger, + + result = RUN_ALL_TESTS(); + return result; +} + + +#endif // MIR_MESH_TEST_H_ diff --git a/src/axom/mir/tests/mir_smoke.cpp b/src/axom/mir/tests/mir_smoke.cpp new file mode 100644 index 0000000000..c71092b54d --- /dev/null +++ b/src/axom/mir/tests/mir_smoke.cpp @@ -0,0 +1,38 @@ +// 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) + +#ifndef MIR_SMOKE_H_ +#define MIR_SMOKE_H_ + +#include "gtest/gtest.h" + +#include "axom/slic.hpp" +#include "axom/mir.hpp" + + +TEST(mir,smoke) +{ + SLIC_INFO("Running smoke test for MIR component"); + + EXPECT_TRUE( true ); + EXPECT_EQ( 0, 0 ); +} + + +//---------------------------------------------------------------------- + +int main(int argc, char* argv[]) +{ + int result = 0; + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::UnitTestLogger logger; // create & initialize test logger, + + result = RUN_ALL_TESTS(); + return result; +} + + +#endif // MIR_SMOKE_H_ diff --git a/src/axom/mir/tests/mir_utilities.cpp b/src/axom/mir/tests/mir_utilities.cpp new file mode 100644 index 0000000000..6d5a6c2dde --- /dev/null +++ b/src/axom/mir/tests/mir_utilities.cpp @@ -0,0 +1,30 @@ +// 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) + +#ifndef MIR_UTILITIES_TEST_H_ +#define MIR_UTILITIES_TEST_H_ + +#include "gtest/gtest.h" + +#include "axom/slic.hpp" +#include "axom/mir.hpp" + + + +//---------------------------------------------------------------------- + +int main(int argc, char* argv[]) +{ + int result = 0; + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::UnitTestLogger logger; // create & initialize test logger, + + result = RUN_ALL_TESTS(); + return result; +} + + +#endif // MIR_UTILITIES_TEST_H_ diff --git a/src/axom/slam/Map.hpp b/src/axom/slam/Map.hpp index 35d9c55090..60da43c560 100644 --- a/src/axom/slam/Map.hpp +++ b/src/axom/slam/Map.hpp @@ -15,7 +15,6 @@ #include #include -#include #include "axom/core.hpp" #include "axom/slic.hpp" @@ -483,23 +482,16 @@ bool Map::isValid(bool verboseOutput) const } - if(verboseOutput) + if(verboseOutput && !bValid) { std::stringstream sstr; sstr << "\n*** Detailed results of isValid on the map.\n"; - if(bValid) - { - sstr << "Map was valid." << std::endl; - } - else - { - sstr << "Map was NOT valid.\n" + sstr << "Map was NOT valid.\n" << sstr.str() << std::endl; - } - std::cout << sstr.str() << std::endl; + SLIC_DEBUG( sstr.str() ); } return bValid; @@ -510,10 +502,11 @@ template void Map::print() const { bool valid = isValid(true); - std::stringstream sstr; if (valid) { + std::stringstream sstr; + if (!m_set) { sstr << "** map is empty."; @@ -534,9 +527,14 @@ void Map::print() const } } } + + SLIC_INFO( sstr.str() ); + } + else + { + SLIC_INFO("Map was not valid."); } - std::cout << sstr.str() << std::endl; } diff --git a/src/docs/dependencies.dot b/src/docs/dependencies.dot index 8e747c824c..5a5b7fef92 100644 --- a/src/docs/dependencies.dot +++ b/src/docs/dependencies.dot @@ -2,6 +2,7 @@ digraph dependencies { quest -> {slam primal mint spin}; {quest slam primal mint spin} -> {slic core}; mint -> sidre [style="dashed"]; + mir -> {slic core slam primal}; spin -> {slam primal}; sidre -> {slic core}; slic -> core; diff --git a/src/docs/doxygen/Doxyfile.in b/src/docs/doxygen/Doxyfile.in index d38a7d5aad..d76fccfee4 100644 --- a/src/docs/doxygen/Doxyfile.in +++ b/src/docs/doxygen/Doxyfile.in @@ -786,6 +786,7 @@ INPUT = @PROJECT_SOURCE_DIR@/axom/mainpage.md \ @PROJECT_SOURCE_DIR@/axom/mint/mesh \ @PROJECT_SOURCE_DIR@/axom/mint/utils \ @PROJECT_SOURCE_DIR@/axom/mint/execution \ + @PROJECT_SOURCE_DIR@/axom/mir/README.md \ @PROJECT_SOURCE_DIR@/axom/primal/README.md \ @PROJECT_SOURCE_DIR@/axom/primal/geometry \ @PROJECT_SOURCE_DIR@/axom/primal/operators \ diff --git a/src/index.rst b/src/index.rst index c4ac47919c..393cba7c82 100644 --- a/src/index.rst +++ b/src/index.rst @@ -61,6 +61,7 @@ for Axom software components: Spin (Spatial indexes) Quest (Querying on surface tool) Mint (Mesh data model) + Mir (Material interface reconstruction) Primal (Computational geometry primitives) -------------------------- @@ -71,6 +72,7 @@ Source Code Documentation * `Core <../../../doxygen/axom_doxygen/html/coretop.html>`_ * `Lumberjack <../../../doxygen/axom_doxygen/html/lumberjacktop.html>`_ * `Mint <../../../doxygen/axom_doxygen/html/minttop.html>`_ + * `Mir <../../../doxygen/axom_doxygen/html/mirtop.html>`_ * `Primal <../../../doxygen/axom_doxygen/html/primaltop.html>`_ * `Quest <../../../doxygen/axom_doxygen/html/questtop.html>`_ * `Sidre <../../../doxygen/axom_doxygen/html/sidretop.html>`_ @@ -86,6 +88,7 @@ Dependencies between modules are as follows: - Slic optionally depends on Lumberjack - Slam, Spin, Primal, Mint, Quest, and Sidre depend on Slic - Mint optionally depends on Sidre +- Mir depends on Slic, Slam and Primal - Quest depends on Slam, Spin, Primal, and Mint The figure below summarizes the dependencies between the modules. Solid links