Merge pull request #1069 from LLNL/feature/kweiss/import-into-samplin…

…g-shaper

Import existing volume fractions into sampling shaper

RELEASE-NOTES.md

@@ -54,6 +54,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
   instance.
 - Multimat: adds an overload of `MultiMat::setCellMatRel()` that supports setting a
   multi-material relation in a compressed sparse-row (CSR) representation.
+- Quest: Adds ability to import volume fractions into `SamplingShaper` before processing `Klee` input
 
 ### Changed
 - Updates blt submodule to HEAD of develop on 24Jan2023

src/axom/quest/IntersectionShaper.hpp

@@ -311,9 +311,7 @@ class IntersectionShaper : public Shaper
                      sidre::MFEMSidreDataCollection* dc)
     : Shaper(shapeSet, dc)
   {
-#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE)
     m_free_mat_name = "free";
-#endif
   }
 
   //@{
@@ -322,6 +320,26 @@ class IntersectionShaper : public Shaper
   void setLevel(int level) { m_level = level; }
 
   void setExecPolicy(int policy) { m_execPolicy = (ExecPolicy)policy; }
+
+  /*!
+   * \brief Set the name of the material used to account for free volume fractions.
+   * \param name The new name of the material. This name cannot contain 
+   *             underscores and it cannot be set once shaping has started.
+   * \note This should not be called once any shaping has occurred.
+   */
+  void setFreeMaterialName(const std::string& name)
+  {
+    if(name.find("_") != std::string::npos)
+    {
+      SLIC_ERROR("The free material name cannot contain underscores.");
+    }
+    if(m_num_elements > 0)
+    {
+      SLIC_ERROR(
+        "The free material name cannot be set once shaping has occurred.");
+    }
+    m_free_mat_name = name;
+  }
   //@}
 
   /*!
@@ -974,27 +992,8 @@ class IntersectionShaper : public Shaper
 
   // Switch back to public. This is done here because the CUDA compiler
   // does not like the following template functions to be private.
-public:
-  /*!
-   * \brief Set the name of the material used to account for free volume fractions.
-   * \param name The new name of the material. This name cannot contain 
-   *             underscores and it cannot be set once shaping has started.
-   * \note This should not be called once any shaping has occurred.
-   */
-  void setFreeMaterialName(const std::string& name)
-  {
-    if(name.find("_") != std::string::npos)
-    {
-      SLIC_ERROR("The free material name cannot contain underscores.");
-    }
-    if(m_num_elements > 0)
-    {
-      SLIC_ERROR(
-        "The free material name cannot be set once shaping has occurred.");
-    }
-    m_free_mat_name = name;
-  }
 
+public:
   /*!
    * \brief Make a new grid function that contains all of the free space not
    *        occupied by existing materials.
@@ -1894,6 +1893,8 @@ class IntersectionShaper : public Shaper
   int m_level {DEFAULT_CIRCLE_REFINEMENT_LEVEL};
   double m_revolvedVolume {DEFAULT_REVOLVED_VOLUME};
   int m_num_elements {0};
+  std::string m_free_mat_name;
+
   double* m_hex_volumes {nullptr};
   double* m_overlap_volumes {nullptr};
 #if defined(AXOM_USE_RAJA) && defined(AXOM_USE_UMPIRE)
@@ -1906,7 +1907,6 @@ class IntersectionShaper : public Shaper
 
   std::vector<mfem::GridFunction*> m_vf_grid_functions;
   std::vector<std::string> m_vf_material_names;
-  std::string m_free_mat_name;
 #endif
 };
 

src/axom/quest/SamplingShaper.hpp

@@ -296,12 +296,12 @@ class SamplingShaper : public Shaper
   void prepareShapeQuery(klee::Dimensions shapeDimension,
                          const klee::Shape& shape) override
   {
-    const auto& shapeName = shape.getName();
+    internal::ScopedLogLevelChanger logLevelChanger(
+      this->isVerbose() ? slic::message::Debug : slic::message::Warning);
 
     SLIC_INFO(axom::fmt::format("{:-^80}", " Generating the octree "));
 
-    internal::ScopedLogLevelChanger logLevelChanger(
-      this->isVerbose() ? slic::message::Debug : slic::message::Warning);
+    const auto& shapeName = shape.getName();
 
     switch(shapeDimension)
     {
@@ -347,13 +347,13 @@ class SamplingShaper : public Shaper
 
   void runShapeQuery(const klee::Shape& shape) override
   {
+    internal::ScopedLogLevelChanger logLevelChanger(
+      this->isVerbose() ? slic::message::Debug : slic::message::Warning);
+
     SLIC_INFO(axom::fmt::format(
       "{:-^80}",
       axom::fmt::format(" Querying the octree for shape '{}'", shape.getName())));
 
-    internal::ScopedLogLevelChanger logLevelChanger(
-      this->isVerbose() ? slic::message::Debug : slic::message::Warning);
-
     switch(getShapeDimension())
     {
     case klee::Dimensions::Two:
@@ -369,15 +369,15 @@ class SamplingShaper : public Shaper
   {
     using axom::utilities::string::rsplitN;
 
+    internal::ScopedLogLevelChanger logLevelChanger(
+      this->isVerbose() ? slic::message::Debug : slic::message::Warning);
+
     const auto& shapeName = shape.getName();
     SLIC_INFO(axom::fmt::format(
       "{:-^80}",
       axom::fmt::format("Applying replacement rules over for shape '{}'",
                         shapeName)));
 
-    internal::ScopedLogLevelChanger logLevelChanger(
-      this->isVerbose() ? slic::message::Debug : slic::message::Warning);
-
     // Get inout qfunc for this shape
     auto* shapeQFunc =
       m_inoutShapeQFuncs.Get(axom::fmt::format("inout_{}", shapeName));
@@ -456,6 +456,60 @@ class SamplingShaper : public Shaper
   //@}
 
 public:
+  /**
+   * \brief Import an initial set of material volume fractions before shaping
+   *
+   * \param [in] initialGridFuncions The input data as a map from material names to grid functions
+   * 
+   * The imported grid functions are interpolated at quadrature points and registered
+   * with the supplied names as material-based quadrature fields
+   */
+  void importInitialVolumeFractions(
+    const std::map<std::string, mfem::GridFunction*>& initialGridFunctions)
+  {
+    internal::ScopedLogLevelChanger logLevelChanger(
+      this->isVerbose() ? slic::message::Debug : slic::message::Warning);
+
+    auto* mesh = m_dc->GetMesh();
+
+    // ensure we have a starting quadrature field for the positions
+    if(!m_inoutShapeQFuncs.Has("positions"))
+    {
+      shaping::generatePositionsQFunction(mesh,
+                                          m_inoutShapeQFuncs,
+                                          m_quadratureOrder);
+    }
+    auto* positionsQSpace = m_inoutShapeQFuncs.Get("positions")->GetSpace();
+
+    // Interpolate grid functions at quadrature points & register material quad functions
+    // assume all elements have same integration rule
+    for(auto& entry : initialGridFunctions)
+    {
+      const auto& name = entry.first;
+      auto* gf = entry.second;
+
+      SLIC_INFO(
+        axom::fmt::format("Importing volume fraction field for '{}' material",
+                          name));
+
+      if(gf == nullptr)
+      {
+        SLIC_WARNING(axom::fmt::format(
+          "Skipping missing volume fraction field for material '{}'",
+          name));
+        continue;
+      }
+
+      auto* matQFunc = new mfem::QuadratureFunction(*positionsQSpace);
+      const auto& ir = matQFunc->GetSpace()->GetIntRule(0);
+      const auto* interp = gf->FESpace()->GetQuadratureInterpolator(ir);
+      interp->Values(*gf, *matQFunc);
+
+      const auto matName = axom::fmt::format("mat_inout_{}", name);
+      m_inoutMaterialQFuncs.Register(matName, matQFunc, true);
+    }
+  }
+
   void adjustVolumeFractions() override
   {
     internal::ScopedLogLevelChanger logLevelChanger(

src/axom/quest/Shaper.cpp

@@ -159,6 +159,9 @@ void Shaper::loadShapeInternal(const klee::Shape& shape,
 {
   using axom::utilities::string::endsWith;
 
+  internal::ScopedLogLevelChanger logLevelChanger(
+    this->isVerbose() ? slic::message::Debug : slic::message::Warning);
+
   SLIC_INFO(axom::fmt::format(
     "{:-^80}",
     axom::fmt::format(" Loading shape '{}' ", shape.getName())));

src/axom/quest/detail/shaping/shaping_helpers.cpp

@@ -127,10 +127,7 @@ void generatePositionsQFunction(mfem::Mesh* mesh,
   inoutQFuncs.Register("positions", pos_coef, true);
 }
 
-/**
- * Compute volume fractions function for shape on a grid of resolution \a gridRes
- * in region defined by bounding box \a queryBounds
- */
+/// Generate a volume fraction from a quadrature field for \a matField using FCT
 void computeVolumeFractions(const std::string& matField,
                             mfem::DataCollection* dc,
                             QFunctionCollection& inoutQFuncs,
@@ -163,24 +160,35 @@ void computeVolumeFractions(const std::string& matField,
                               dim,
                               NE));
 
-  // Project QField onto volume fractions field
+  // Access or create a registered volume fraction grid function from the data collection
+  mfem::FiniteElementSpace* fes = nullptr;
+  mfem::GridFunction* volFrac = nullptr;
+  if(dc->HasField(volFracName))
+  {
+    volFrac = dc->GetField(volFracName);
+    fes = volFrac->FESpace();
+  }
+  else
+  {
+    const auto basis = mfem::BasisType::Positive;
+    auto* fec = new mfem::L2_FECollection(outputOrder, dim, basis);
+    fes = new mfem::FiniteElementSpace(mesh, fec);
+    volFrac = new mfem::GridFunction(fes);
+    volFrac->MakeOwner(fec);
 
-  mfem::L2_FECollection* fec =
-    new mfem::L2_FECollection(outputOrder, dim, mfem::BasisType::Positive);
-  mfem::FiniteElementSpace* fes = new mfem::FiniteElementSpace(mesh, fec);
-  mfem::GridFunction* volFrac = new mfem::GridFunction(fes);
-  volFrac->MakeOwner(fec);
-  dc->RegisterField(volFracName, volFrac);
+    dc->RegisterField(volFracName, volFrac);
+  }
 
+  // Project QField onto volume fractions field using flux corrected transport (FCT)
+  // to keep the range of values between 0 and 1
   axom::utilities::Timer timer(true);
   {
-    mfem::MassIntegrator mass_integrator;  // use the default for exact integration; lower for approximate
-
+    mfem::MassIntegrator mass_integrator;
     mfem::QuadratureFunctionCoefficient qfc(*inout);
     mfem::DomainLFIntegrator rhs(qfc);
 
-    // assume all elts are the same
-    const auto& ir = inout->GetSpace()->GetIntRule(0);
+    const auto& ir =
+      inout->GetSpace()->GetIntRule(0);  // assume all elements are the same
     rhs.SetIntRule(&ir);
 
     mfem::DenseMatrix m;
@@ -200,9 +208,6 @@ void computeVolumeFractions(const std::string& matField,
       rhs.AssembleRHSElementVect(*el, *T, b);
       mInv.Factor(m);
 
-      // Use FCT limiting -- similar to Remap
-      // Limit the function to be between 0 and 1
-      // Q: Is there a better way limiting algorithm for this?
       if(one.Size() != b.Size())
       {
         one.SetSize(b.Size());

src/axom/quest/detail/shaping/shaping_helpers.hpp

@@ -55,8 +55,16 @@ void generatePositionsQFunction(mfem::Mesh* mesh,
                                 int sampleRes);
 
 /**
- * Compute volume fractions function for shape on a grid of resolution \a gridRes
- * in region defined by bounding box \a queryBounds
+ * \brief Compute volume fractions for a given material using its associated quadrature function
+ *
+ * \param [in] matField The name of the material
+ * \param [in] dc The DataCollection containing the specified material
+ * \param [in] inoutQFuncs A collection of quadrature functions containing the quadrature
+ * values associated with the specified material
+ * \param [in] outputOrder The order the grid function that we're generating
+ *
+ * The generated grid function will be prefixed by `vol_frac_`
+ * 
  */
 void computeVolumeFractions(const std::string& matField,
                             mfem::DataCollection* dc,