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             : #pragma once
      11             : 
      12             : #include "MooseUtils.h"
      13             : #include "MooseHashing.h"
      14             : #include "libmesh/point.h"
      15             : 
      16             : /// An unordered map indexed by Point, eg 3 floating point numbers
      17             : /// Because floating point rounding errors can affect the hashing, eg the binning of values,
      18             : /// we may need to look near a Point/key for where the initial value was stored
      19             : struct PointIndexedMap
      20             : {
      21      116681 :   PointIndexedMap(const Point & mesh_max_coords)
      22      116681 :   {
      23      466724 :     for (auto i : make_range(LIBMESH_DIM))
      24      350043 :       normalization(i) =
      25      700086 :           MooseUtils::absoluteFuzzyEqual(mesh_max_coords(i), 0) ? 1 : mesh_max_coords(i);
      26      116681 :   }
      27             : 
      28             :   /**
      29             :    * Normalize, expand then round a point to create local bins
      30             :    * - normalize the point by dividing by the meshes max (non-zero) absolute coordinates
      31             :    * - multiply the result by 1e12 and round, which keeps a fixed number of digits
      32             :    */
      33     7781120 :   inline Point getRoundedPoint(const Point & p) const
      34             :   {
      35             :     // We cant support coordinates too far from the normalization point
      36             :     mooseAssert(p(0) / normalization(0) < 1 && p(1) / normalization(1) < 1 &&
      37             :                     p(2) / normalization(2) < 1,
      38             :                 "Point coordinates for indexing must be normalized below 1. "
      39             :                 "Point: "
      40             :                     << p << " normalization: " << normalization);
      41             : 
      42     7781120 :     std::vector<Real> rounded_p{round(p(0) / normalization(0) * pow(10, 12)),
      43     7781120 :                                 round(p(1) / normalization(1) * pow(10, 12)),
      44    15562240 :                                 round(p(2) / normalization(2) * pow(10, 12))};
      45    15562240 :     return Point(rounded_p[0], rounded_p[1], rounded_p[2]);
      46     7781120 :   }
      47             : 
      48             :   // This needs to be moved into libMesh according to Fande
      49             :   // The hash helps sort the points in buckets
      50             :   // Hashing the floats also leads to floating point comparison errors in their own way
      51             :   // We hash a rounding of the float instead
      52             :   struct hash_point
      53             :   {
      54     8274018 :     std::size_t operator()(const Point & p) const
      55             :     {
      56     8274018 :       std::size_t seed = 0;
      57     8274018 :       Moose::hash_combine(seed, p(0), p(1), p(2));
      58     8274018 :       return seed;
      59             :     }
      60             :   };
      61             : 
      62             :   /// The container indexed by points
      63             :   std::unordered_map<Point, Number, hash_point> base_map;
      64             : 
      65             :   /// Normalization factors used to scale points before insertions/comparisons
      66             :   Point normalization;
      67             : 
      68      462974 :   Number & operator[](Point p)
      69             :   {
      70      462974 :     auto it = find(p);
      71      462974 :     if (it != end())
      72      136689 :       return const_cast<Number &>(it->second);
      73             :     else
      74             :     {
      75      326285 :       const auto rp = getRoundedPoint(p);
      76      326285 :       return base_map[rp];
      77             :     }
      78             :   }
      79             : 
      80      360676 :   unsigned int hasKey(Point p)
      81             :   {
      82      360676 :     auto it = find(p);
      83      360676 :     if (it != end())
      84       95476 :       return true;
      85             :     else
      86      265200 :       return false;
      87             :   }
      88             : 
      89     1089471 :   std::unordered_map<Point, Number, hash_point>::const_iterator find(const Point & pt) const
      90             :   {
      91     1089471 :     const auto rp = getRoundedPoint(pt);
      92             : 
      93     1089471 :     auto it = base_map.find(rp);
      94     1089471 :     if (it != base_map.end())
      95      293793 :       return it;
      96             :     else
      97             :     {
      98             :       // Search the point with coordinates rounded down to origin
      99      795678 :       auto rounded = base_map.find(rp);
     100      795678 :       if (rounded != base_map.end())
     101           0 :         return rounded;
     102             : 
     103             :       // To store the output of tentative searches
     104      795678 :       std::unordered_map<Point, Number, hash_point>::const_iterator out;
     105             : 
     106             :       // Search in all directions around
     107     2386974 :       for (int i : make_range(2))
     108     4774008 :         for (int j : make_range(2))
     109     9548016 :           for (int k : make_range(2))
     110     6365364 :             if (attempt_find(pt, 2 * i - 1, 2 * j - 1, 2 * k - 1, out))
     111          60 :               return out;
     112             : 
     113      795618 :       return base_map.end();
     114             :     }
     115             :   }
     116             : 
     117     6365364 :   bool attempt_find(const Point & p,
     118             :                     Real dx,
     119             :                     Real dy,
     120             :                     Real dz,
     121             :                     std::unordered_map<Point, Number, hash_point>::const_iterator & out) const
     122             :   {
     123             :     // Some epsilon to move the point by.
     124             :     // There is no multiplicative epsilon that will reliably change the last decimal
     125     6365364 :     const auto eps = 1e-13;
     126    12730728 :     const auto shifted_rp = getRoundedPoint(p + Point(dx * eps * normalization(0),
     127     6365364 :                                                       dy * eps * normalization(1),
     128    12730728 :                                                       dz * eps * normalization(2)));
     129     6365364 :     out = base_map.find(shifted_rp);
     130     6365364 :     if (out != base_map.end())
     131          60 :       return true;
     132     6365304 :     return false;
     133             :   }
     134             : 
     135     1089471 :   std::unordered_map<Point, Number, hash_point>::const_iterator end() const
     136             :   {
     137     1089471 :     return base_map.end();
     138             :   }
     139             : };