Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #ifdef MOOSE_MFEM_ENABLED
      11             : 
      12             : #include "MFEMQuadratureFunctionCoefficientBase.h"
      13             : #include "MooseError.h"
      14             : 
      15             : #include "libmesh/ignore_warnings.h"
      16             : #include "mfem.hpp"
      17             : #include "libmesh/restore_warnings.h"
      18             : 
      19             : namespace
      20             : {
      21             : 
      22             : bool
      23      554840 : SamePoint(const mfem::IntegrationPoint & a, const mfem::IntegrationPoint & b)
      24             : {
      25      554840 :   constexpr mfem::real_t tol = 1e-12;
      26     2219342 :   auto close = [](mfem::real_t x, mfem::real_t y) { return !(std::abs(x - y) > tol); };
      27      554840 :   return close(a.x, b.x) && close(a.y, b.y) && close(a.z, b.z) && close(a.weight, b.weight);
      28             : }
      29             : }
      30             : 
      31             : void
      32      554834 : MFEMQuadratureFunctionCoefficientBase::CheckIntegrationRule(const mfem::QuadratureFunction & qf,
      33             :                                                             mfem::ElementTransformation & T,
      34             :                                                             const mfem::IntegrationPoint & ip) const
      35             : {
      36      554834 :   const mfem::QuadratureSpaceBase & qspace = *qf.GetSpace();
      37      554834 :   const int el_idx = qspace.GetEntityIndex(T);
      38             :   // Entity not in this space (e.g. a boundary evaluation); mfem's Eval returns zero, nothing to
      39             :   // check.
      40      554834 :   if (el_idx < 0)
      41           0 :     return;
      42             : 
      43      554834 :   const mfem::IntegrationRule & stored_rule = qspace.GetIntRule(el_idx);
      44             :   // Fast path: the consuming integrator's point at this index coincides with the stored rule's,
      45             :   // so the rules match.
      46      554834 :   if (ip.index < stored_rule.Size() && SamePoint(stored_rule.IntPoint(ip.index), ip))
      47      554832 :     return;
      48             : 
      49             :   // Rules differ. Recover the quadrature order the coefficient should have used by finding the
      50             :   // lowest order whose rule for this geometry reproduces the consuming integrator's point at
      51             :   // ip.index. This runs only on the error path.
      52           2 :   const int stored_order = qspace.GetOrder();
      53           2 :   const mfem::Geometry::Type geom = qspace.GetGeometry(el_idx);
      54           2 :   int suggested_order = -1;
      55             :   // The consuming integrator's order is not known here, but is realistically bounded. The upper
      56             :   // limit is just a generous finite cap so this error-path search always terminates: a multiple
      57             :   // of the stored order, plus a constant so low-order stored rules still search far enough. If the
      58             :   // true order somehow exceeds the cap, suggested_order stays -1 and a generic message is emitted.
      59           6 :   for (int order = 0; order <= 2 * stored_order + 64; ++order)
      60             :   {
      61           6 :     const mfem::IntegrationRule & candidate = mfem::IntRules.Get(geom, order);
      62           6 :     if (ip.index < candidate.Size() && SamePoint(candidate.IntPoint(ip.index), ip))
      63             :     {
      64           2 :       suggested_order = order;
      65           2 :       break;
      66             :     }
      67             :   }
      68             : 
      69           2 :   if (suggested_order >= 0)
      70           2 :     mooseError("MFEM quadrature function '",
      71           2 :                _name,
      72             :                "' stores values on the order-",
      73             :                stored_order,
      74             :                " quadrature rule (",
      75           2 :                stored_rule.Size(),
      76             :                " points on this element), but it is being evaluated by an integrator using a "
      77             :                "different quadrature rule (",
      78           2 :                mfem::IntRules.Get(geom, suggested_order).Size(),
      79             :                " points). The stored values are indexed by quadrature point, so the orders must "
      80             :                "match. Set 'order = ",
      81             :                suggested_order,
      82             :                "' on '",
      83           2 :                _name,
      84             :                "'.");
      85             :   else
      86           0 :     mooseError("MFEM quadrature function '",
      87           0 :                _name,
      88             :                "' stores values on the order-",
      89             :                stored_order,
      90             :                " quadrature rule, but it is being evaluated by an integrator using a different "
      91             :                "quadrature rule. The stored values are indexed by quadrature point, so the orders "
      92             :                "must match; adjust 'order' on '",
      93           0 :                _name,
      94             :                "' to match the consuming integrator.");
      95             : }
      96             : 
      97             : #endif