Line data    Source code

       1             : /**********************************************************************/
       2             : /*                     DO NOT MODIFY THIS HEADER                      */
       3             : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
       4             : /*                                                                    */
       5             : /*            Copyright 2017 Battelle Energy Alliance, LLC            */
       6             : /*                        ALL RIGHTS RESERVED                         */
       7             : /**********************************************************************/
       8             : 
       9             : #include "PointLocatorRegularOrthogonalData.h"
      10             : #include "MooseError.h"
      11             : #include "libmesh/mesh_base.h"
      12             : #include "libmesh/elem.h"
      13             : 
      14          95 : PointLocatorRegularOrthogonalData::PointLocatorRegularOrthogonalData(
      15             :     const std::vector<unsigned int> & cell_count,
      16             :     const Point & min_corner,
      17             :     const Point & max_corner,
      18          95 :     const MeshBase & mesh)
      19          95 :   : _dim(cell_count.size()),
      20          95 :     _cell_count(cell_count),
      21          95 :     _min_corner(min_corner),
      22             :     _cell_size(max_corner - min_corner)
      23             : {
      24             :   // determine total number of cells and their size
      25             :   unsigned int num_root_elems = 1;
      26         352 :   for (unsigned int i = 0; i < _dim; ++i)
      27             :   {
      28         257 :     if (_cell_count[i] == 0)
      29           0 :       mooseError("Cannot have zero cells in any spatial direction");
      30             : 
      31         257 :     num_root_elems *= _cell_count[i];
      32         257 :     _cell_size(i) /= _cell_count[i];
      33             :   }
      34             : 
      35             :   // initialize cell storage
      36          95 :   _root_elems.assign(num_root_elems, nullptr);
      37             : 
      38             :   // sort all level 0 elements into the table
      39          95 :   auto el = mesh.level_elements_begin(0);
      40          95 :   const auto end_el = mesh.level_elements_end(0);
      41     1176861 :   for (; el != end_el; ++el)
      42             :   {
      43      588383 :     const Elem * elem = *el;
      44      588383 :     const unsigned int index = rootElementIndex(elem->vertex_average());
      45      588383 :     _root_elems[index] = elem;
      46             :   }
      47             : 
      48             :   // sanity check
      49      588478 :   for (auto el : _root_elems)
      50      588383 :     if (el == nullptr)
      51           0 :       mooseError("Found a null root element");
      52          95 : }
      53             : 
      54             : const Elem *
      55    73438975 : PointLocatorRegularOrthogonalData::rootElement(const Point & p, Point & el_pos) const
      56             : {
      57             :   unsigned int index = 0;
      58             :   const Point p0 = p - _min_corner;
      59             : 
      60   288605561 :   for (unsigned int i = 0; i < _dim; ++i)
      61             :   {
      62   215191744 :     const int n = p0(i) / _cell_size(i);
      63             : 
      64             :     // outside of mesh, return null
      65   215191744 :     if (n < 0 || n >= _cell_count[i])
      66             :       return nullptr;
      67             : 
      68             :     // cell internal coordinates from 0..1 used for bisecting refined cells
      69   215166586 :     el_pos(i) = p0(i) / _cell_size(i) - n;
      70             : 
      71             :     // construct root element index
      72   215166586 :     index = index * _cell_count[i] + n;
      73             :   }
      74             : 
      75    73413817 :   return _root_elems[index];
      76             : }
      77             : 
      78             : unsigned int
      79      588383 : PointLocatorRegularOrthogonalData::rootElementIndex(const Point & p) const
      80             : {
      81             :   unsigned int index = 0;
      82     2310876 :   for (unsigned int i = 0; i < _dim; ++i)
      83             :   {
      84     1722493 :     const int n = (p(i) - _min_corner(i)) / _cell_size(i);
      85             : 
      86     1722493 :     if (n < 0 || n >= _cell_count[i])
      87           0 :       mooseError("Point not inside regular orthogonal mesh");
      88             : 
      89     1722493 :     index = index * _cell_count[i] + n;
      90             :   }
      91             : 
      92      588383 :   return index;
      93             : }