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, const int dim)
      24             : {
      25      554840 :   constexpr mfem::real_t tol = 1e-12;
      26             :   mfem::real_t pa[3], pb[3];
      27      554840 :   a.Get(pa, dim);
      28      554840 :   b.Get(pb, dim);
      29     1775100 :   for (int d = 0; d < dim; ++d)
      30     1220266 :     if (std::abs(pa[d] - pb[d]) > tol)
      31           6 :       return false;
      32      554834 :   return std::abs(a.weight - b.weight) < tol;
      33             : }
      34             : }
      35             : 
      36             : void
      37      554834 : MFEMQuadratureFunctionCoefficientBase::CheckIntegrationRule(const mfem::QuadratureFunction & qf,
      38             :                                                             mfem::ElementTransformation & T,
      39             :                                                             const mfem::IntegrationPoint & ip) const
      40             : {
      41      554834 :   const mfem::QuadratureSpaceBase & qspace = *qf.GetSpace();
      42      554834 :   const int el_idx = qspace.GetEntityIndex(T);
      43             :   // Entity not in this space (e.g. a boundary evaluation)
      44      554834 :   if (el_idx < 0)
      45           0 :     return;
      46             : 
      47      554834 :   const mfem::Geometry::Type geom = qspace.GetGeometry(el_idx);
      48      554834 :   const int dim = mfem::Geometry::Dimension[geom];
      49      554834 :   const mfem::IntegrationRule & stored_rule = qspace.GetIntRule(el_idx);
      50             :   // Fast path: the consuming integrator's point at this index coincides with the stored rule's,
      51             :   // so the rules match.
      52      554834 :   if (ip.index < stored_rule.Size() && SamePoint(stored_rule.IntPoint(ip.index), ip, dim))
      53      554832 :     return;
      54             : 
      55             :   // Rules differ. Recover the quadrature order the coefficient should have used by finding the
      56             :   // lowest order whose rule for this geometry reproduces the consuming integrator's point at
      57             :   // ip.index. This runs only on the error path.
      58           2 :   const int stored_order = qspace.GetOrder();
      59           2 :   int suggested_order = -1;
      60             :   // The consuming integrator's order is not known here, but is realistically bounded. The upper
      61             :   // limit is just a generous finite cap so this error-path search always terminates: a multiple
      62             :   // of the stored order, plus a constant so low-order stored rules still search far enough. If the
      63             :   // true order somehow exceeds the cap, suggested_order stays -1 and a generic message is emitted.
      64           6 :   for (int order = 0; order <= 2 * stored_order + 64; ++order)
      65             :   {
      66           6 :     const mfem::IntegrationRule & candidate = mfem::IntRules.Get(geom, order);
      67           6 :     if (ip.index < candidate.Size() && SamePoint(candidate.IntPoint(ip.index), ip, dim))
      68             :     {
      69           2 :       suggested_order = order;
      70           2 :       break;
      71             :     }
      72             :   }
      73             : 
      74           2 :   if (suggested_order >= 0)
      75           2 :     mooseError("MFEM quadrature function '",
      76           2 :                _name,
      77             :                "' stores values on the order-",
      78             :                stored_order,
      79             :                " quadrature rule (",
      80           2 :                stored_rule.Size(),
      81             :                " points on this element), but it is being evaluated by an integrator using a "
      82             :                "different quadrature rule (",
      83           2 :                mfem::IntRules.Get(geom, suggested_order).Size(),
      84             :                " points). The stored values are indexed by quadrature point, so the orders must "
      85             :                "match. Set 'order = ",
      86             :                suggested_order,
      87             :                "' on '",
      88           2 :                _name,
      89             :                "'.");
      90             :   else
      91           0 :     mooseError("MFEM quadrature function '",
      92           0 :                _name,
      93             :                "' stores values on the order-",
      94             :                stored_order,
      95             :                " quadrature rule, but it is being evaluated by an integrator using a different "
      96             :                "quadrature rule. The stored values are indexed by quadrature point, so the orders "
      97             :                "must match; adjust 'order' on '",
      98           0 :                _name,
      99             :                "' to match the consuming integrator.");
     100             : }
     101             : 
     102             : #endif