bb_grid_map.cpp

Go to the documentation of this file.

// Copyright (c) 2010-2026, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
 
/* The BBoxTensorGridMap class is adapted from similar functionality in the
gslib library. Below is the gslib license and copyright statement:
 
Copyright (c) 2008-2024, UCHICAGO ARGONNE, LLC.
 
The UChicago Argonne, LLC as Operator of Argonne National
Laboratory holds copyright in the Software. The copyright holder
reserves all rights except those expressly granted to licensees,
and U.S. Government license rights.
 
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
 
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the disclaimer below.
 
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the disclaimer (as noted below)
in the documentation and/or other materials provided with the
distribution.
 
3. Neither the name of ANL nor the names of its contributors
may be used to endorse or promote products derived from this software
without specific prior written permission.
 
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
UCHICAGO ARGONNE, LLC, THE U.S. DEPARTMENT OF
ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
 
#include "bb_grid_map.hpp"
 
#include <limits>
#include <cmath>
#include <algorithm>
 
namespace mfem
{
 
using namespace std;
 
BBoxTensorGridMap::BBoxTensorGridMap(Mesh &mesh, int nx)
{
   GridFunction *nodes = mesh.GetNodes();
   const int nel = mesh.GetNE();
   sdim = mesh.SpaceDimension();
   Vector elmin(nel*sdim), elmax(nel*sdim);
   elmin = numeric_limits<real_t>::max();
   elmax = -numeric_limits<real_t>::max();
   if (!nodes)
   {
      Array<int> verts;
      real_t *coord;
      // create bounding boxes from vertex coordinates
      for (int e = 0; e < nel; e++)
      {
         mesh.GetElementVertices(e, verts);
         for (int v = 0; v < verts.Size(); v++)
         {
            coord = mesh.GetVertex(verts[v]);
            for (int d = 0; d < sdim; d++)
            {
               elmin(d*nel + e) = min(elmin(d*nel + e), coord[d]);
               elmax(d*nel + e) = max(elmax(d*nel + e), coord[d]);
            }
         }
      }
   }
   else
   {
      int nref = 3;
      nodes->GetElementBounds(elmin, elmax, nref);
   }
   Array<int> nx_arr(sdim);
   nx_arr = nx;
   Setup(elmin, elmax, nel, nx_arr, false);
}
 
BBoxTensorGridMap::BBoxTensorGridMap(Vector &elmin,
                                     Vector &elmax,
                                     int nel,
                                     int sdim_,
                                     int n,
                                     bool by_max_size)
{
   sdim = sdim_;
   MFEM_VERIFY(0 < sdim && sdim <= 3,
               "BBoxTensorGridMap only supports spatial dimensions 1, 2, and 3.");
   if (nel > 0)
   {
      MFEM_VERIFY(elmin.Size() == sdim * nel && elmax.Size() == sdim * nel,
                  "Element bounds size must match dim * nel.");
   }
   Array<int> nx_arr(sdim);
   nx_arr = n;
   Setup(elmin, elmax, nel, nx_arr, by_max_size);
}
 
BBoxTensorGridMap::BBoxTensorGridMap(Vector &elmin, Vector &elmax,
                                     int nel, int sdim_,
                                     Array<int> &nx,
                                     bool by_max_size)
{
   sdim = sdim_;
   Setup(elmin, elmax, nel, nx, by_max_size);
}
 
void BBoxTensorGridMap::Setup(Vector &elmin, Vector &elmax,
                              int nel, Array<int> &nx, bool by_max_size)
{
   MFEM_VERIFY(0 < sdim && sdim <= 3,
               "BBoxTensorGridMap only supports spatial dimensions 1, 2, and 3.");
   MFEM_VERIFY(nx.Size() == sdim,
               "BBoxTensorGridMap requires nx to have the same size as the number of dimensions.");
   if (nel > 0)
   {
      MFEM_VERIFY(elmin.Size() == sdim * nel && elmax.Size() == sdim * nel,
                  "Element bounds size must match dim * nel.");
   }
   lmap_bnd_min.SetSize(sdim);
   lmap_bnd_max.SetSize(sdim);
   lmap_fac.SetSize(sdim);
   lmap_nx.SetSize(sdim);
   lmap_nx = nx;
 
   if (by_max_size)
   {
      MFEM_VERIFY(nx[0] >= 0,
                  "BBoxTensorGridMap requires a nonnegative max-size hint.");
   }
   else
   {
      for (int d = 0; d < nx.Size(); d++)
      {
         MFEM_VERIFY(nx[d] > 0,
                     "BBoxTensorGridMap requires positive number of divisions in each dimension.");
      }
   }
   if (nel == 0)
   {
      lmap_bnd_min = 0.0;
      lmap_bnd_max = 1.0;
      if (by_max_size) { lmap_nx = 1; }
      SetGridFac(lmap_fac, lmap_nx, lmap_bnd_min, lmap_bnd_max);
 
      lmap_nxd = lmap_nx[0];
      for (int d = 1; d < sdim; d++)
      {
         lmap_nxd *= lmap_nx[d];
      }
 
      lgrid_map.SetSize(lmap_nxd + 1);
      lgrid_map = lmap_nxd + 1;
      return;
   }
 
   for (int d = 0; d < sdim; d++)
   {
      Vector elmind(elmin.GetData() + d*nel, nel);
      Vector elmaxd(elmax.GetData() + d*nel, nel);
      lmap_bnd_min[d] = elmind.Min();
      lmap_bnd_max[d] = elmaxd.Max();
   }
 
   Array<int> elmin_h, elmax_h;
   unsigned int store_size;
   if (by_max_size)
   {
      int nmax = nx[0];
      int nlow = 1, nhigh = nmax > nel ? ceil(pow(nmax - nel, 1.0 / sdim)) : 1;
      int size_low = 2 + nel;
      int size = 0;
      while (nhigh - nlow > 1)
      {
         int nmid = nlow + (nhigh - nlow) / 2;
         int nmd = nmid;
         for (int d = 1; d < sdim; d++)
         {
            nmd *= nmid;
         }
         lmap_nx = nmid;
         SetGridFac(lmap_fac, lmap_nx, lmap_bnd_min, lmap_bnd_max);
         size = nmd + 1 + GetGridCountAndRange(lmap_nx, lmap_fac,
                                               lmap_bnd_min, lmap_bnd_max,
                                               elmin, elmax,
                                               elmin_h, elmax_h);
         if (size <= nmax) { nlow = nmid; size_low = size; }
         else { nhigh = nmid; }
      }
      lmap_nx = nlow;
      lmap_nxd = nlow;
      for (int d = 1; d < sdim; d++)
      {
         lmap_nxd *= nlow;
      }
      store_size = size_low;
      SetGridFac(lmap_fac, lmap_nx, lmap_bnd_min, lmap_bnd_max);
      if (size != size_low)
      {
         GetGridCountAndRange(lmap_nx, lmap_fac,
                              lmap_bnd_min, lmap_bnd_max,
                              elmin, elmax,
                              elmin_h, elmax_h);
      }
   }
   else
   {
      SetGridFac(lmap_fac, lmap_nx, lmap_bnd_min, lmap_bnd_max);
 
      lmap_nxd = lmap_nx[0];
      for (int d = 1; d < sdim; d++)
      {
         lmap_nxd *= lmap_nx[d];
      }
 
      // Grid cell ranges for each element in each direction
      store_size = lmap_nxd + 1 + GetGridCountAndRange(lmap_nx, lmap_fac,
                                                       lmap_bnd_min,
                                                       lmap_bnd_max,
                                                       elmin, elmax,
                                                       elmin_h, elmax_h);
   }
 
   lgrid_map.SetSize(store_size);
   lgrid_map[0] = lmap_nxd + 1;
 
   Array<unsigned int> grid_el_count(lmap_nxd);
   grid_el_count = 0;
 
   for (int e = 0; e < nel; e++)
   {
      int klim = sdim < 3 ? 1 : (elmax_h[2*nel+e]-elmin_h[2*nel+e]);
      int jlim = sdim < 2 ? 1 : (elmax_h[1*nel+e]-elmin_h[1*nel+e]);
      int ilim = (elmax_h[0*nel+e]-elmin_h[0*nel+e]);
      for (int k = 0; k < klim; k++)
      {
         int koff = sdim < 3 ? 0 :
                    (elmin_h[2*nel + e] + k) * lmap_nx[0] * lmap_nx[1];
         for (int j = 0; j < jlim; j++)
         {
            int joff = sdim < 2 ? 0 : (elmin_h[1*nel + e] + j) * lmap_nx[0];
            for (int i = 0; i < ilim; i++)
            {
               int ioff = elmin_h[e] + i;
               int idx = ioff + joff + koff;
               grid_el_count[idx]++;
            }
         }
      }
   }
 
   for (unsigned int e = 0; e < lmap_nxd; e++)
   {
      lgrid_map[e + 1] = lgrid_map[e] + grid_el_count[e];
   }
 
   for (int e = 0; e < nel; e++)
   {
      int klim = sdim < 3 ? 1 : (elmax_h[2*nel+e]-elmin_h[2*nel+e]);
      int jlim = sdim < 2 ? 1 : (elmax_h[1*nel+e]-elmin_h[1*nel+e]);
      int ilim = (elmax_h[0*nel+e]-elmin_h[0*nel+e]);
      for (int k = 0; k < klim; k++)
      {
         int koff = sdim < 3 ? 0 :
                    (elmin_h[2*nel+e] + k) * lmap_nx[0] * lmap_nx[1];
         for (int j = 0; j < jlim; j++)
         {
            int joff = sdim < 2 ? 0 : (elmin_h[1*nel + e] + j) * lmap_nx[0];
            for (int i = 0; i < ilim; i++)
            {
               int ioff = elmin_h[e] + i;
               int idx = ioff + joff + koff;
               lgrid_map[lgrid_map[idx+1]-grid_el_count[idx]]=e;
               grid_el_count[idx]--;
            }
         }
      }
   }
}
 
Array<int> BBoxTensorGridMap::GridCellToElements(int i) const
{
   MFEM_ASSERT(i >= 0 && (unsigned int)i < lmap_nxd,
               "Access element " << i << " of local grid with cells = "
               << lmap_nxd);
   int start = lgrid_map[i];
   int end = lgrid_map[i + 1];
   Array<int> elements(end - start);
   for (int j = start; j < end; j++)
   {
      elements[j - start] = lgrid_map[j];
   }
   return elements;
}
 
int BBoxTensorGridMap::GetGridCellFromPoint(Vector &xyz) const
{
   MFEM_ASSERT(xyz.Size() == sdim,
               "Point must have the same dimension as the grid.");
   int sum = 0;
   for (int d = sdim-1; d >= 0; --d)
   {
      if (xyz(d) < lmap_bnd_min(d) || xyz(d) > lmap_bnd_max(d))
      {
         return -1; // Point is outside the bounds of the grid
      }
      sum *= lmap_nx[d];
      int i = (int)floor((xyz(d) - lmap_bnd_min(d)) * lmap_fac[d]);
      sum += i < 0 ? 0 : (lmap_nx[d] - 1 < i ? lmap_nx[d] - 1 : i);
   }
   return sum;
}
 
Array<int> BBoxTensorGridMap::MapPointToElements(Vector &xyz) const
{
   MFEM_ASSERT(xyz.Size() == sdim,
               "Point must have the same dimension as the grid.");
   int cell = GetGridCellFromPoint(xyz);
   if (cell < 0)
   {
      return Array<int>(); // Point is outside the bounds of the tensor grid
   }
   return GridCellToElements(cell);
}
 
void BBoxTensorGridMap::GetGridRange(const int d, const Array<int> &lh_n,
                                     const Vector &lh_fac,
                                     const Vector &lh_bnd_min,
                                     const real_t &xmin, const real_t &xmax,
                                     int &imin, int &imax)
{
   // Use a half-open interval [imin, imax) for the covered grid-cell range.
   // If xmin is exactly on a grid boundary, use the cell on the right/high
   // side. If xmax is exactly on a grid boundary, stop before the cell on the
   // right/high side.
   int i0 = floor( (xmin - lh_bnd_min[d]) * lh_fac[d] );
   int i1 = ceil ( (xmax - lh_bnd_min[d]) * lh_fac[d] );
   imin = std::clamp(i0, 0, lh_n[d] - 1);
   imax = std::clamp(i1, imin + 1, lh_n[d]);
}
 
void BBoxTensorGridMap::SetGridFac(Vector &lh_fac, const Array<int> &nx,
                                   const Vector &lh_bnd_min,
                                   const Vector &lh_bnd_max)
{
   int dim = lh_bnd_min.Size();
   for (int d = 0; d < dim; d++)
   {
      real_t length = lh_bnd_max[d] - lh_bnd_min[d];
      if (length > 0.0)
      {
         lh_fac[d] = nx[d] / length;
      }
      else
      {
         lh_fac[d] = 0.0;
      }
   }
}
 
int BBoxTensorGridMap::GetGridCountAndRange(const Array<int> &lh_n,
                                            const Vector &lh_fac,
                                            const Vector &lh_bnd_min,
                                            const Vector &lh_bnd_max,
                                            const Vector &elmin,
                                            const Vector &elmax,
                                            Array<int> &elmin_h,
                                            Array<int> &elmax_h)
{
   int count = 0;
   const int dim = lh_bnd_min.Size();
   const int nel = elmin.Size()/dim;
   elmin_h.SetSize(dim * nel);
   elmax_h.SetSize(dim * nel);
   for (int i = 0; i < nel; i++)
   {
      int count_el = 1;
      for (int d = 0; d < dim; d++)
      {
         GetGridRange(d, lh_n, lh_fac, lh_bnd_min,
                      elmin[d*nel + i], elmax[d*nel + i],
                      elmin_h[d*nel + i], elmax_h[d*nel + i]);
         int imax = elmax_h[d*nel + i];
         int imin = elmin_h[d*nel + i];
         count_el *= (imax - imin);
      }
      count += count_el;
   }
   return count;
}
 
} // namespace mfem