gslib.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.
 
#include "gslib.hpp"
#include "geom.hpp"
#include "../mesh/bb_grid_map.hpp"
#include "../general/forall.hpp"
#include <algorithm> // std::fill
 
#ifdef MFEM_USE_GSLIB
 
/* The class GlobalBBoxTensorGridMap, and certain methods including
dbl_range_diag_expand_2, dbl_range_diag_expand_3, bbox_2_tfm, bbox_3_tfm,
dbl_range_expand, and obbox{surf/edge}_calc_{D}, were adapted from 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.
*/
 
// Ignore warnings from the gslib header (GCC version)
#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-function"
#endif
 
#define CODE_INTERNAL 0
#define CODE_BORDER 1
#define CODE_NOT_FOUND 2
 
// External GSLIB header (the MFEM header is gslib.hpp)
namespace gslib
{
#include "gslib.h"
#ifndef GSLIB_RELEASE_VERSION //gslib v1.0.7
#define GSLIB_RELEASE_VERSION 10007
#endif
static_assert(std::is_same_v<uint,unsigned int>,
              "GSLIB's integer-type, 'uint', defined in gslib.h, must be the same as 'unsigned int'!");
 
extern "C" {
   struct hash_data_3
   {
      ulong hash_n;
      struct dbl_range bnd[3];
      double fac[3];
      uint *offset;
   };
 
   struct hash_data_2
   {
      ulong hash_n;
      struct dbl_range bnd[2];
      double fac[2];
      uint *offset;
   };
 
   struct findpts_dummy_ms_data
   {
      unsigned int *nsid;
      double       *distfint;
   };
 
   struct findpts_data_3
   {
      struct crystal cr;
      struct findpts_local_data_3 local;
      struct hash_data_3 hash;
      struct array savpt;
      struct findpts_dummy_ms_data fdms;
      uint   fevsetup;
   };
 
   struct findpts_data_2
   {
      struct crystal cr;
      struct findpts_local_data_2 local;
      struct hash_data_2 hash;
      struct array savpt;
      struct findpts_dummy_ms_data fdms;
      uint   fevsetup;
   };
} //extern C
 
} //namespace gslib
 
#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
 
namespace mfem
{
 
// Static helper to map device-side find-points outputs to MFEM reference-space
// data. nvcc does not allow lambdas in non-public members.
static void FindPointsDeviceSetCode(int pts_cnt,
                                    int ddim,
                                    real_t rbtol,
                                    real_t bdr_t,
                                    unsigned int *d_gsl_code,
                                    const real_t *d_gsl_ref,
                                    const real_t *d_gsl_dist,
                                    const unsigned int *d_gsl_elem,
                                    real_t *d_gsl_mfem_ref,
                                    unsigned int *d_gsl_mfem_elem)
{
   mfem::forall(pts_cnt, [=] MFEM_HOST_DEVICE (int index)
   {
      if (d_gsl_code[index] == CODE_NOT_FOUND)
      {
         return;
      }
      d_gsl_mfem_elem[index] = d_gsl_elem[index];
 
      bool internal = true;
      for (int d = 0; d < ddim; d++)
      {
         const real_t r_val = d_gsl_ref[index * ddim + d];
         const real_t val = 0.5 * (r_val + 1.0);
         d_gsl_mfem_ref[index * ddim + d] = val;
         if (val < rbtol || val > 1.0 - rbtol)
         {
            internal = false;
         }
      }
 
      const int setcode = internal ? CODE_INTERNAL : CODE_BORDER;
      d_gsl_code[index] = (setcode == CODE_BORDER &&
                           d_gsl_dist[index] > bdr_t)
                          ? CODE_NOT_FOUND : setcode;
   });
}
 
// Static helper: scatter interpolated values back into field_out.
// nvcc does not allow lambdas in non-public members.
static void InterpolateDeviceScatter(int nlocal,
                                     const int *d_index_temp,
                                     const real_t *d_interp_vals,
                                     real_t *d_field_out,
                                     int interp_offset,
                                     int ncomp,
                                     int pts_cnt,
                                     int ordering)
{
   if (nlocal == 0 || ncomp == 0) { return; }
 
   mfem::forall(nlocal*ncomp, [=] MFEM_HOST_DEVICE (int k)
   {
      const int j = k % nlocal;
      const int i = k / nlocal;
      const int pt_index = d_index_temp[j];
      const int idx = (ordering == Ordering::byNODES) ?
                      pt_index + i*pts_cnt :
                      pt_index*ncomp + i;
      d_field_out[idx] = d_interp_vals[j + interp_offset*i];
   });
}
 
FindPointsGSLIB::FindPointsGSLIB()
   : mesh(NULL),
     fec_map_lin(NULL),
     fdataD(NULL), cr(NULL), gsl_comm(NULL),
     dim(-1), spacedim(-1), points_cnt(-1), setupflag(false),
     default_interp_value(0),
     avgtype(AvgType::ARITHMETIC), bdr_tol(1e-8)
{
   mesh_split.SetSize(4);
   ir_split.SetSize(4);
   ir_split_sol.SetSize(4);
   fes_rst_map.SetSize(4);
   gf_rst_map.SetSize(4);
   for (int i = 0; i < mesh_split.Size(); i++)
   {
      mesh_split[i] = nullptr;
      ir_split[i] = nullptr;
      ir_split_sol[i] = nullptr;
      fes_rst_map[i] = nullptr;
      gf_rst_map[i] = nullptr;
   }
 
   gsl_comm = new gslib::comm;
#ifdef MFEM_USE_MPI
   if (!Mpi::IsInitialized()) { Mpi::Init(); }
   MPI_Comm comm = MPI_COMM_WORLD;
   comm_init(gsl_comm, comm);
#else
   comm_init(gsl_comm, 0);
#endif
}
 
FindPointsGSLIB::FindPointsGSLIB(Mesh &mesh_in,
                                 const double bbox_rel_size_inc,
                                 const double newt_tol, const int npt_max)
   : FindPointsGSLIB()
{
   Setup(mesh_in, bbox_rel_size_inc, newt_tol, npt_max);
}
 
FindPointsGSLIB::~FindPointsGSLIB()
{
   FreeData();
#ifdef MFEM_USE_MPI
   if (!Mpi::IsFinalized())  // currently segfaults inside gslib otherwise
#endif
   {
      comm_free(gsl_comm);
      delete gsl_comm;
   }
   for (int i = 0; i < mesh_split.Size(); i++)
   {
      if (mesh_split[i]) { delete mesh_split[i]; mesh_split[i] = nullptr; }
      if (ir_split[i]) { delete ir_split[i]; ir_split[i] = nullptr; }
      if (ir_split_sol[i]) { delete ir_split_sol[i]; ir_split_sol[i] = nullptr; }
      if (fes_rst_map[i]) { delete fes_rst_map[i]; fes_rst_map[i] = nullptr; }
      if (gf_rst_map[i]) { delete gf_rst_map[i]; gf_rst_map[i] = nullptr; }
   }
   if (fec_map_lin) { delete fec_map_lin; fec_map_lin = nullptr; }
}
 
#ifdef MFEM_USE_MPI
FindPointsGSLIB::FindPointsGSLIB(MPI_Comm comm_)
   : mesh(NULL),
     fec_map_lin(NULL),
     fdataD(NULL), cr(NULL), gsl_comm(NULL),
     dim(-1), spacedim(-1), points_cnt(-1), setupflag(false),
     default_interp_value(0),
     avgtype(AvgType::ARITHMETIC), bdr_tol(1e-8)
{
   mesh_split.SetSize(4);
   ir_split.SetSize(4);
   ir_split_sol.SetSize(4);
   fes_rst_map.SetSize(4);
   gf_rst_map.SetSize(4);
   for (int i = 0; i < mesh_split.Size(); i++)
   {
      mesh_split[i] = nullptr;
      ir_split[i] = nullptr;
      ir_split_sol[i] = nullptr;
      fes_rst_map[i] = nullptr;
      gf_rst_map[i] = nullptr;
   }
 
   gsl_comm = new gslib::comm;
   comm_init(gsl_comm, comm_);
}
 
FindPointsGSLIB::FindPointsGSLIB(ParMesh &mesh_in,
                                 const double bbox_rel_size_inc,
                                 const double newt_tol, const int npt_max)
   : FindPointsGSLIB(mesh_in.GetComm())
{
   Setup(mesh_in, bbox_rel_size_inc, newt_tol, npt_max);
}
#endif
 
void FindPointsGSLIB::SetupCrystal()
{
   if (cr == NULL)
   {
      cr = new gslib::crystal;
      crystal_init(cr, gsl_comm);
   }
}
 
void FindPointsGSLIB::FreeCrystal()
{
   if (cr == NULL) { return; }
 
#ifdef MFEM_USE_MPI
   if (!Mpi::IsFinalized())  // currently segfaults inside gslib otherwise
#endif
   {
      crystal_free(cr);
   }
   delete cr;
   cr = NULL;
   DEV.cr = NULL;
}
 
void FindPointsGSLIB::Setup(Mesh &m, const double bbox_rel_size_inc,
                            const double newt_tol,
                            const int npt_max)
{
   MFEM_VERIFY(m.GetNodes() != nullptr, "Mesh nodes are required.");
   const int meshOrder = m.GetNodes()->FESpace()->GetMaxElementOrder();
 
   // call FreeData if FindPointsGSLIB::Setup has been called already
   if (m.Dimension() != m.SpaceDimension())
   {
      SetupSurf(m, bbox_rel_size_inc, newt_tol);
      return;
   }
   if (setupflag) { FreeData(); }
 
   SetupCrystal();
 
   mesh = &m;
   dim  = mesh->Dimension();
   spacedim = dim;
   const unsigned dof1D = meshOrder + 1;
 
   SetupSplitMeshesAndIntegrationRules(meshOrder);
 
   GetNodalValues(mesh->GetNodes(), gsl_mesh);
 
   mesh_points_cnt = gsl_mesh.Size()/dim;
 
   DEV.local_hash_size = mesh_points_cnt;
   DEV.dof1d = (int)dof1D;
   if (dim == 2)
   {
      unsigned nr[2] = { dof1D, dof1D };
      unsigned mr[2] = { 2*dof1D, 2*dof1D };
      double * const elx[2] =
      {
         mesh_points_cnt == 0 ? nullptr : &gsl_mesh(0),
         mesh_points_cnt == 0 ? nullptr : &gsl_mesh(mesh_points_cnt)
      };
      fdataD = findpts_setup_2(gsl_comm, elx, nr, NE_split_total, mr,
                               bbox_rel_size_inc, DEV.local_hash_size,
                               mesh_points_cnt, npt_max, newt_tol);
   }
   else
   {
      unsigned nr[3] = { dof1D, dof1D, dof1D };
      unsigned mr[3] = { 2*dof1D, 2*dof1D, 2*dof1D };
      double * const elx[3] =
      {
         mesh_points_cnt == 0 ? nullptr : &gsl_mesh(0),
         mesh_points_cnt == 0 ? nullptr : &gsl_mesh(mesh_points_cnt),
         mesh_points_cnt == 0 ? nullptr : &gsl_mesh(2*mesh_points_cnt)
      };
      fdataD = findpts_setup_3(gsl_comm, elx, nr, NE_split_total, mr,
                               bbox_rel_size_inc, DEV.local_hash_size,
                               mesh_points_cnt, npt_max, newt_tol);
   }
   setupflag = true;
}
 
 
/* Calculates the diagonal length of the bounding box and expands its bounds by
 * 0.5*len*tol at both its min and max values.
 * Returns the length of the diagonal (could be used for expanding obboxes).
 */
double dbl_range_diag_expand_2(struct gslib::dbl_range *b, double tol)
{
   double l[2] = { b[0].max-b[0].min, b[1].max-b[1].min };
   double len = sqrt(l[0]*l[0] + l[1]*l[1])*0.5*tol;
   for (int i=0; i<2; i++)
   {
      b[i].min = b[i].min - len;
      b[i].max = b[i].max + len;
   }
   return len;
}
 
double dbl_range_diag_expand_3(struct gslib::dbl_range *b, double tol)
{
   double l[3] = { b[0].max-b[0].min, b[1].max-b[1].min, b[2].max-b[2].min };
   double len = sqrt(l[0]*l[0] + l[1]*l[1] + l[2]*l[2])*0.5*tol;
   for (int i=0; i<3; i++)
   {
      b[i].min = b[i].min - len;
      b[i].max = b[i].max + len;
   }
   return len;
}
 
static void bbox_2_tfm(double *tfm, const double x0[2], const double Ji[4],
                       const double *x, const double *y, unsigned n)
{
   unsigned i;
   for (i=0; i<n; ++i)
   {
      const double dx = x[i]-x0[0], dy = y[i]-x0[1];
      tfm[  i] = Ji[0]*dx + Ji[1]*dy;
      tfm[n+i] = Ji[2]*dx + Ji[3]*dy;
   }
}
 
static void bbox_3_tfm(double *tfm, const double x0[3], const double Ji[9],
                       const double *x, const double *y, const double *z,
                       unsigned n)
{
   unsigned i;
   for (i=0; i<n; ++i)
   {
      const double dx = x[i]-x0[0], dy = y[i]-x0[1], dz = z[i]-x0[2];
      tfm[    i] = Ji[0]*dx + Ji[1]*dy + Ji[2]*dz;
      tfm[  n+i] = Ji[3]*dx + Ji[4]*dy + Ji[5]*dz;
      tfm[2*n+i] = Ji[6]*dx + Ji[7]*dy + Ji[8]*dz;
   }
}
 
static struct gslib::dbl_range dbl_range_expand(struct gslib::dbl_range b,
                                                double tol)
{
   double a = (b.min+b.max)/2, l = (b.max-b.min)*(1+tol)/2;
   struct gslib::dbl_range m;
   m.min = a-l, m.max = a+l;
   return m;
}
 
static MFEM_HOST_DEVICE inline void MapSplitTriangleQuadToTriangle(
   const int tri_id, const double u, const double v, double &tx, double &ty)
{
   const double N0 = (1.0-u)*(1.0-v);
   const double N1 = u*(1.0-v);
   const double N2 = u*v;
   const double N3 = (1.0-u)*v;
 
   // Must match the quad -> triangle split in SetupSplitMeshes.
   const double vx[3][4] =
   {
      {0.0, 0.5,       1.0/3.0, 0.0    },
      {0.5, 1.0,       0.5,     1.0/3.0},
      {0.0, 1.0/3.0,   0.5,     0.0    }
   };
   const double vy[3][4] =
   {
      {0.0, 0.0,       1.0/3.0, 0.5    },
      {0.0, 0.0,       0.5,     1.0/3.0},
      {0.5, 1.0/3.0,   0.5,     1.0    }
   };
 
   tx = N0*vx[tri_id][0] + N1*vx[tri_id][1]
        + N2*vx[tri_id][2] + N3*vx[tri_id][3];
   ty = N0*vy[tri_id][0] + N1*vy[tri_id][1]
        + N2*vy[tri_id][2] + N3*vy[tri_id][3];
}
 
static void VerifyAABBPadLayout(const Vector *aabb_sz_inc, const uint nel,
                                const int sd)
{
   if (!aabb_sz_inc) { return; }
 
   const int sz = aabb_sz_inc->Size();
   MFEM_VERIFY(sz == 1 || sz == (int)nel || sz == sd || sz == (int)nel*sd,
               "Invalid aabb_sz_inc length for SetupSurfWithAABBExpansion: "
               "expected 1, NE, SpaceDim, or NE*SpaceDim.");
}
 
static double GetAABBPad(const Vector *aabb_sz_inc, const int aabb_sz_inc_size,
                         const uint nel, const int sd, const uint e,
                         const int d)
{
   if (!aabb_sz_inc) { return 0.0; }
 
   MFEM_ASSERT(aabb_sz_inc_size == 1 || aabb_sz_inc_size == (int)nel ||
               aabb_sz_inc_size == sd || aabb_sz_inc_size == (int)nel*sd,
               "Invalid aabb_sz_inc layout.");
 
   double s = 0.0;
   if (aabb_sz_inc_size == 1) { s = (*aabb_sz_inc)(0); }
   else if (aabb_sz_inc_size == (int)nel) { s = (*aabb_sz_inc)((int)e); }
   else if (aabb_sz_inc_size == sd) { s = (*aabb_sz_inc)(d); }
   else { s = (*aabb_sz_inc)((int)e*sd + d); }
 
   MFEM_VERIFY(s >= 0.0,
               "aabb_sz_inc absolute AABB expansion must be non-negative.");
   return 0.5*s;
}
 
void obboxsurf_calc_3(Vector &bb,
                      const double *const elx[3],
                      const unsigned n,
                      uint nel,
                      const unsigned m,
                      const double tol,
                      const bool store_obb)
{
   auto h_bb = bb.HostWrite();
   const double *x = elx[0], *y = elx[1], *z = elx[2];
   const int n_el_ents = store_obb ? 18 : 6;
   // 3(c0) + 3(aabb_min) + 3(aabb_max) + 9(A)
 
   const unsigned n2 = n*n;
   const unsigned lbsize = gslib::lob_bnd_size(n, m);
 
   // 2*n for tensor_ig2
   // 2*m*(n+m+1) for lob_bnd_2 (always greater than 2*n for n and m >=1)
   // gslib::gll_lag_size(n) for lag_setup and then lag
   unsigned wsize = std::max(2*m*(n+m+1), store_obb?gslib::gll_lag_size(n):0u);
 
   Vector datavec(lbsize + wsize + (store_obb ? 2*n : 0));
   double *data = datavec.GetData();
 
   double *const I0 = store_obb ? data : nullptr; // basis and derivative at r=0
   double *const lob_bnd_data = store_obb ? (data + 2*n) : data;
   double *const work = lob_bnd_data + lbsize;
 
   gslib::lob_bnd_setup(lob_bnd_data, n, m); // setup machinery for bounding
 
   if (!store_obb)
   {
      for (int ie = 0; ie < (int)nel; ie++, x += n2, y += n2, z += n2)
      {
         struct gslib::dbl_range ab[3];
         ab[0] = gslib::lob_bnd_2(lob_bnd_data,n,m, lob_bnd_data,n,m, x, work);
         ab[1] = gslib::lob_bnd_2(lob_bnd_data,n,m, lob_bnd_data,n,m, y, work);
         ab[2] = gslib::lob_bnd_2(lob_bnd_data,n,m, lob_bnd_data,n,m, z, work);
 
         dbl_range_diag_expand_3(ab, tol);
 
         h_bb[n_el_ents*ie + 0] = ab[0].min;
         h_bb[n_el_ents*ie + 1] = ab[1].min;
         h_bb[n_el_ents*ie + 2] = ab[2].min;
         h_bb[n_el_ents*ie + 3] = ab[0].max;
         h_bb[n_el_ents*ie + 4] = ab[1].max;
         h_bb[n_el_ents*ie + 5] = ab[2].max;
      }
      return;
   }
 
   gslib::lagrange_fun *const lag = gslib::gll_lag_setup(work, n);
   lag(I0, work, n, 1, 0);
 
   for (int ie = 0; (unsigned)ie < nel; ie++,x+=n2,y+=n2,z+=n2)
   {
      struct gslib::dbl_range ab[3];
      struct gslib::dbl_range tb[3];
      double x0[3], tv[9], A[9];
 
      /*
         * Find the center of the element (r=0 ref. coord.) in physical space
         * and store in x0.
         * tv[0], tv[1], tv[2]: kept empty at this point for convenience.
         * tv[3], tv[4]: dx/dr, dx/ds
         * tv[5], tv[6]: dy/dr, dy/ds
         * tv[7], tv[8]: dz/dr, dz/ds
         */
      x0[0] = gslib::tensor_ig2(tv+3, I0,n, I0,n, x, work);
      x0[1] = gslib::tensor_ig2(tv+5, I0,n, I0,n, y, work);
      x0[2] = gslib::tensor_ig2(tv+7, I0,n, I0,n, z, work);
 
      // tangent vector 1 moved to tv[0], tv[1], tv[2]
      tv[0] = tv[3], tv[1] = tv[5], tv[2] = tv[7];
      // tangent vector 2 moved to tv[3], tv[4], tv[5]
      tv[3] = tv[4], tv[4] = tv[6], tv[5] = tv[8];
      // normal vector to the plane formed by t1 and t2 (cross product)
      // is stored in tv[6], tv[7], tv[8]
      tv[6] = tv[1]*tv[5] - tv[2]*tv[4];
      tv[7] = tv[2]*tv[3] - tv[0]*tv[5];
      tv[8] = tv[0]*tv[4] - tv[1]*tv[3];
      // normalize the normal vector
      const double nmag  = sqrt(tv[6]*tv[6] + tv[7]*tv[7] + tv[8]*tv[8]);
      tv[6] = tv[6]/nmag;
      tv[7] = tv[7]/nmag;
      tv[8] = tv[8]/nmag;
 
      // Rodrigues formula to compute the rotation matrix
      // Axis of rotation is n x [0,0,1] = [n_2, -n_1, 0], and we must
      // normalize it
      double nmag2 = tv[6]*tv[6] + tv[7]*tv[7];
      if (nmag2 > 0)
      {
         nmag2 = sqrt(nmag2);
         tv[7] = tv[7]/nmag2;
         tv[6] = tv[6]/nmag2;
      }
      double kx = tv[7];
      double ky = -tv[6];
      double kz = 0.0;
 
      double ct = tv[8];
      double st = nmag2; //1.0 - ct*ct;
 
      // row-major
      A[0] = 1.0 + st*0.0 + (1.0-ct)*(-ky*ky-kz*kz);
      A[1] = 0.0 + st*(0.0) + (1.0-ct)*(kx*ky);
      A[2] = 0.0 + st*(ky) + (1.0-ct)*(kx*kz);
 
      A[3] = 0.0 + st*(0.0) + (1.0-ct)*(kx*ky);
      A[4] = 1.0 + st*(0.0) + (1.0-ct)*(-kx*kx-kz*kz);
      A[5] = 0.0 + st*(-kx) + (1.0-ct)*(ky*kz);
 
      A[6] = 0.0 + st*(-ky) + (1.0-ct)*(kx*kz);
      A[7] = 0.0 + st*(kx) + (1.0-ct)*(ky*kz);
      A[8] = 1.0 + st*(0.0) + (1.0-ct)*(-kx*kx-ky*ky);
 
      // compute bounds in physical space
      ab[0] = gslib::lob_bnd_2(lob_bnd_data,n,m, lob_bnd_data,n,m, x, work);
      ab[1] = gslib::lob_bnd_2(lob_bnd_data,n,m, lob_bnd_data,n,m, y, work);
      ab[2] = gslib::lob_bnd_2(lob_bnd_data,n,m, lob_bnd_data,n,m, z, work);
      // expand bounding box based on (tol*diagonal_length) in each direction
      // to avoid 0 extent in 1 direction.
      double aabb_diag_len = dbl_range_diag_expand_3(ab, tol);
      h_bb[n_el_ents*ie + 3] = ab[0].min;
      h_bb[n_el_ents*ie + 4] = ab[1].min;
      h_bb[n_el_ents*ie + 5] = ab[2].min;
      h_bb[n_el_ents*ie + 6] = ab[0].max;
      h_bb[n_el_ents*ie + 7] = ab[1].max;
      h_bb[n_el_ents*ie + 8] = ab[2].max;
 
      // rotate to align normal with z-axis
      Array<double> xtfm(3*n2);
      bbox_3_tfm(xtfm.GetData(), x0,A, x,y,z,n2);
      // The rotated z-coords are used to calculate z-bounds.
      tb[2] = gslib::lob_bnd_2(lob_bnd_data,n,m, lob_bnd_data,n,m,
                               xtfm.GetData()+2*n2, work);
 
      tb[2].min -= aabb_diag_len;
      tb[2].max += aabb_diag_len;
 
      // Also apply A to the tangent vectors, which allows us to
      // calculate the Jacobian matrix at the rotated element center.
      // NOTE that the z components of the rotated tangent vectors will
      // become zero, since the normal vector is aligned with the z-axis.
      double J[4], Ji[4];
      J[0] = A[0]*tv[0] + A[1]*tv[1] + A[2]*tv[2]; // rotated dx/dr
      J[1] = A[0]*tv[3] + A[1]*tv[4] + A[2]*tv[5]; // rotated dx/ds
      J[2] = A[3]*tv[0] + A[4]*tv[1] + A[5]*tv[2]; // rotated dy/dr
      J[3] = A[3]*tv[3] + A[4]*tv[4] + A[5]*tv[5]; // rotated dy/ds
      DenseMatrix JM(J, 2, 2);
      DenseMatrix JiM(Ji, 2, 2);
      CalcInverse(JM, JiM);
 
      // Now transform the already rotated x,y coordinates according to Ji to
      // their reference space.
      // Important to note that the nodes used here already have the element
      // center at (0,0). Hence, Ji can be directly applied to them.
      for (unsigned i=0; i<n2; ++i)
      {
         const double xt = xtfm[i], yt = xtfm[n2+i];
         xtfm[   i] = Ji[0]*xt + Ji[1]*yt;
         xtfm[n2+i] = Ji[2]*xt + Ji[3]*yt;
      }
      // Bound these reference space xy coordinates
      tb[0] = gslib::lob_bnd_2(lob_bnd_data,n,m, lob_bnd_data,n,m,
                               xtfm.GetData(), work);
      tb[1] = gslib::lob_bnd_2(lob_bnd_data,n,m, lob_bnd_data,n,m,
                               xtfm.GetData()+n2, work);
      // Expand the bounds based on the tol
      tb[0] = dbl_range_expand(tb[0],tol);
      tb[1] = dbl_range_expand(tb[1],tol);
 
      /* We now have a BB whose bounds represent bounds of a OBB around the
         * original element.
         *
         * We calculate the center of the OBB in physical space by calculating
         * the center of this BB, which is the same as the displacement needed
         * to move from the element center in the transformed space to the BB center. This displacement is then untransformed by applying (Ji.A)^-1
         * to it, and added to known physical element center.
         *
         * This BB does not necessarily have known fixed size like [-1,1].
         * So, we premultiply a length scaling matrix, say L, to Ji.A to
         * L.Ji.A. This is the total transformation needed to move a physical
         * location that is inside the physical OBB to a location within [-1,1]^3. Any transformed point not in [-1,1]^3 is outside the OBB.
         *
         * It must be noted: this transformation matrix is only applied to
         * points that have been translated by the physical OBB center.
         */
      {
         // The center of the BB in the transformed space
         const double av0 = (tb[0].min+tb[0].max)/2,
                      av1 = (tb[1].min+tb[1].max)/2,
                      av2 = (tb[2].min+tb[2].max)/2;
         // First untransform the x,y coordinates by J to obtain all
         // components in the rotated space
         const double Jav0 = J[0]*av0 + J[1]*av1,
                      Jav1 = J[2]*av0 + J[3]*av1;
         // The physical displacement needed to move from the element center
         // to the OBB center is calculated by "un"rotating {Jav0,Jav1,av2}
         // by applying inverse of A.
         // The physical untransformed OBB center can then be obtained.
 
         h_bb[n_el_ents*ie + 0] = x0[0] + A[0]*Jav0 + A[3]*Jav1 + A[6]*av2;
         h_bb[n_el_ents*ie + 1] = x0[1] + A[1]*Jav0 + A[4]*Jav1 + A[7]*av2;
         h_bb[n_el_ents*ie + 2] = x0[2] + A[2]*Jav0 + A[5]*Jav1 + A[8]*av2;
      }
 
      // Finally, obtain (L.Ji.A) and store it in out->A
      {
         // The scaling matrix L's diagonal terms, needed to scale the
         // transformation to [-1,1]^3.
         const double di0 = 2/(tb[0].max-tb[0].min),
                      di1 = 2/(tb[1].max-tb[1].min),
                      di2 = 2/(tb[2].max-tb[2].min);
 
         // We finally construct the final transformation matrix A=L.Ji.A.
         // This maps a position relative to OBB center to a position in
         // [-1,1]^3, if the position is inside the OBB.
         h_bb[n_el_ents*ie + 9 ]=di0*(Ji[0]*A[0] + Ji[1]*A[3]);
         h_bb[n_el_ents*ie + 10]=di0*(Ji[0]*A[1] + Ji[1]*A[4]);
         h_bb[n_el_ents*ie + 11]=di0*(Ji[0]*A[2] + Ji[1]*A[5]);
         h_bb[n_el_ents*ie + 12]=di1*(Ji[2]*A[0] + Ji[3]*A[3]);
         h_bb[n_el_ents*ie + 13]=di1*(Ji[2]*A[1] + Ji[3]*A[4]);
         h_bb[n_el_ents*ie + 14]=di1*(Ji[2]*A[2] + Ji[3]*A[5]);
         h_bb[n_el_ents*ie + 15]=di2*A[6];
         h_bb[n_el_ents*ie + 16]=di2*A[7];
         h_bb[n_el_ents*ie + 17]=di2*A[8];
      }
   }
}
 
void obboxedge_calc_2(Vector &bb,
                      const double *const elx[2],
                      const unsigned nr,
                      uint nel,
                      const unsigned mr,
                      const double tol,
                      const bool store_obb)
{
   auto h_bb = bb.HostWrite();
   const double *x   = elx[0];
   const double *y   = elx[1];
 
   const int n_el_ents = store_obb ? 10 : 4; // store only AABB if !store_obb
 
   const unsigned lbsize0 = gslib::lob_bnd_size(nr,mr);
   unsigned wsize = std::max(2*nr+2*mr,store_obb ? gslib::gll_lag_size(nr) : 0);
   Vector datavec(lbsize0 + wsize + (store_obb ? 2*nr : 0));
   double *data = datavec.GetData();
 
   double *const I0r = store_obb ? data : nullptr;
   double *const lob_bnd_data_r = store_obb ? data + 2*nr : data;
   double *const work = lob_bnd_data_r + lbsize0;
 
   gslib::lob_bnd_setup(lob_bnd_data_r, nr, mr);
 
   if (!store_obb)
   {
      for (int ie = 0; ie < (int)nel; ie++, x+=nr, y+=nr)
      {
         struct gslib::dbl_range ab[2];
         ab[0] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, x, work);
         ab[1] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, y, work);
 
         dbl_range_diag_expand_2(ab, tol);
 
         h_bb[n_el_ents*ie + 0] = ab[0].min;
         h_bb[n_el_ents*ie + 1] = ab[1].min;
         h_bb[n_el_ents*ie + 2] = ab[0].max;
         h_bb[n_el_ents*ie + 3] = ab[1].max;
      }
      return;
   }
 
   gslib::lagrange_fun *const lag = gslib::gll_lag_setup(work, nr);
   lag(I0r, work, nr,1, 0);
 
   for (int ie = 0; (unsigned)ie < nel; ie++,x+=nr,y+=nr)
   {
      double x0[2], A[4];
      struct gslib::dbl_range ab[2], tb[2];
 
      // Find the physical-space coordinates of center of the element (r=0).
      // A holds the Jacobian/tangent dx/dr and dy/dr
      x0[0] = gslib::tensor_ig1(A,I0r,nr,x);
      x0[1] = gslib::tensor_ig1(A+1,I0r,nr,y);
      // Normalize the tangent and construct the 2x2 rotation matrix A that
      // aligns the tangent with the x-axis. Rows of A are [tangent, normal].
      A[2] = sqrt(A[0]*A[0] + A[1]*A[1]);
      A[0] = A[0]/A[2];
      A[1] = A[1]/A[2];
      A[2] = -A[1];
      A[3] =  A[0];
 
      // Compute AABB in x and y.
      ab[0] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, x, work);
      ab[1] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, y, work);
      // Transform nodes to the reference frame and compute bounds.
      bbox_2_tfm(work, x0, A, x, y, nr);
      tb[0] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, work,    work+2*nr);
      tb[1] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, work+nr, work+2*nr);
 
      // expand AABB
      double aabb_diag_len = dbl_range_diag_expand_2(ab, tol);
 
      // The center of the OBB in the rotated frame.
      const double av0 = (tb[0].min+tb[0].max)/2,
                   av1 = (tb[1].min+tb[1].max)/2;
      // A is orthogonal, so A^-1 = A^T. Untransform [av0,av1] back to
      // physical space to get the physical OBB center.
      const double dx0 =  A[0]*av0 - A[1]*av1,
                   dx1 = -A[2]*av0 + A[3]*av1;
      h_bb[n_el_ents*ie + 0] = x0[0] + dx0;
      h_bb[n_el_ents*ie + 1] = x0[1] + dx1;
      h_bb[n_el_ents*ie + 2] = ab[0].min;
      h_bb[n_el_ents*ie + 3] = ab[1].min;
      h_bb[n_el_ents*ie + 4] = ab[0].max;
      h_bb[n_el_ents*ie + 5] = ab[1].max;
 
      // Expand the rotated-frame bounds by aabb_diag_len, then premultiply
      // A by the scaling matrix L = diag(di0,di1) to obtain L*A. This maps
      // a physical displacement from the OBB center to [-1,1]^2.
      tb[0].min -= aabb_diag_len;
      tb[0].max += aabb_diag_len;
      tb[1].min -= aabb_diag_len;
      tb[1].max += aabb_diag_len;
      const double di0 = 2/(tb[0].max-tb[0].min),
                   di1 = 2/(tb[1].max-tb[1].min);
      h_bb[n_el_ents*ie + 6]=di0*A[0];
      h_bb[n_el_ents*ie + 7]=di0*A[1];
      h_bb[n_el_ents*ie + 8]=di1*A[2];
      h_bb[n_el_ents*ie + 9]=di1*A[3];
   }
}
 
void obboxedge_calc_3(Vector &bb,
                      const double *const elx[3],
                      const unsigned nr,
                      uint nel,
                      const unsigned mr,
                      const double tol,
                      const bool store_obb)
{
   auto h_bb = bb.HostWrite();
   const double *x   = elx[0];
   const double *y   = elx[1];
   const double *z   = elx[2];
   const int n_el_ents = store_obb ? 18 : 6;
 
   const unsigned lbsize0 = gslib::lob_bnd_size(nr,mr);
   // 2*mr for lob_bnd_1
   // +3*nr for storing rotated coordinates computed by bbox_3_tfm
   const unsigned wsize = store_obb ?
                          std::max(3*nr+2*mr, gslib::gll_lag_size(nr)) :
                          (2*nr + 2*mr);
 
   Vector datavec(lbsize0 + wsize + (store_obb ? 2*nr : 0));
   double *data = datavec.GetData();
 
   double *const I0r = store_obb ? data : nullptr;
   double *const lob_bnd_data_r = store_obb ? data + 2*nr : data;
   double *const work = lob_bnd_data_r + lbsize0;
 
   gslib::lob_bnd_setup(lob_bnd_data_r, nr, mr);
 
   if (!store_obb)
   {
      for (int ie = 0; ie < (int)nel; ie++, x += nr, y += nr, z += nr)
      {
         struct gslib::dbl_range ab[3];
         ab[0] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, x, work);
         ab[1] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, y, work);
         ab[2] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, z, work);
 
         dbl_range_diag_expand_3(ab, tol);
 
         h_bb[n_el_ents*ie + 0] = ab[0].min;
         h_bb[n_el_ents*ie + 1] = ab[1].min;
         h_bb[n_el_ents*ie + 2] = ab[2].min;
         h_bb[n_el_ents*ie + 3] = ab[0].max;
         h_bb[n_el_ents*ie + 4] = ab[1].max;
         h_bb[n_el_ents*ie + 5] = ab[2].max;
      }
      return;
   }
 
   gslib::lagrange_fun *const lag = gslib::gll_lag_setup(work, nr);
   lag(I0r, work, nr, 1, 0);
 
   for (int ie = 0; (unsigned)ie < nel; ie++,x+=nr,y+=nr,z+=nr)
   {
      double x0[3], A[9], Ai[9];
      struct gslib::dbl_range ab[3], tb[3];
 
      // Find the physical-space coordinates of center of the element (r=0).
      // A holds the Jacobian/tangent dx/dr, dy/dr, and dz/dr.
      x0[0] = gslib::tensor_ig1(A,I0r,nr,x);
      x0[1] = gslib::tensor_ig1(A+1,I0r,nr,y);
      x0[2] = gslib::tensor_ig1(A+2,I0r,nr,z);
 
      // Normalize the tangent vector.
      double nmag  = A[0]*A[0] + A[1]*A[1] + A[2]*A[2];
      if (nmag > 0)
      {
         nmag = sqrt(nmag);
         A[0] = A[0]/nmag;
         A[1] = A[1]/nmag;
         A[2] = A[2]/nmag;
      }
 
      double nmag2 = A[0]*A[0] + A[1]*A[1];
      if (nmag2 > 0)
      {
         nmag2 = sqrt(nmag2);
         A[1] = A[1]/nmag2;
         A[0] = A[0]/nmag2;
      }
      // Rodrigues formula: axis of rotation is tangent x [1,0,0],
      // normalized to k = [kx,ky,kz]. ct and st are cos/sin of the angle.
      double kx = A[1];
      double ky = -A[0];
      double kz = 0.0;
 
      double ct = A[2];
      double st = nmag2; //1.0 - ct*ct;
 
      // Construct the 3x3 rotation matrix A that aligns the tangent with
      // the x-axis (row-major).
      A[0] = 1.0 + st*0.0 + (1.0-ct)*(-ky*ky-kz*kz);
      A[1] = 0.0 + st*(0.0) + (1.0-ct)*(kx*ky);
      A[2] = 0.0 + st*(ky) + (1.0-ct)*(kx*kz);
 
      A[3] = 0.0 + st*(0.0) + (1.0-ct)*(kx*ky);
      A[4] = 1.0 + st*(0.0) + (1.0-ct)*(-kx*kx-kz*kz);
      A[5] = 0.0 + st*(-kx) + (1.0-ct)*(ky*kz);
 
      A[6] = 0.0 + st*(-ky) + (1.0-ct)*(kx*kz);
      A[7] = 0.0 + st*(kx) + (1.0-ct)*(ky*kz);
      A[8] = 1.0 + st*(0.0) + (1.0-ct)*(-kx*kx-ky*ky);
 
      DenseMatrix AM(A, 3, 3);
      DenseMatrix AiM(Ai, 3, 3);
      CalcInverse(AM, AiM);
 
      // Compute AABB in x, y, and z.
      ab[0] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, x, work);
      ab[1] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, y, work);
      ab[2] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, z, work);
      // Transform nodes to the rotated frame and compute bounds.
      bbox_3_tfm(work, x0, A, x, y, z, nr);
      tb[0] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, work,      work+3*nr);
      tb[1] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, work+nr,   work+3*nr);
      tb[2] = gslib::lob_bnd_1(lob_bnd_data_r, nr, mr, work+2*nr, work+3*nr);
 
      // expand AABB
      double aabb_diag_len = dbl_range_diag_expand_3(ab, tol);
 
      // The center of the OBB in the rotated frame.
      const double av0 = (tb[0].min+tb[0].max)/2,
                   av1 = (tb[1].min+tb[1].max)/2,
                   av2 = (tb[2].min+tb[2].max)/2;
      // Untransform [av0,av1,av2] back to physical space using A^-1
      // to obtain the physical OBB center.
      h_bb[n_el_ents*ie + 0] = x0[0] + Ai[0]*av0 + Ai[1]*av1 + Ai[2]*av2;
      h_bb[n_el_ents*ie + 1] = x0[1] + Ai[3]*av0 + Ai[4]*av1 + Ai[5]*av2;
      h_bb[n_el_ents*ie + 2] = x0[2] + Ai[6]*av0 + Ai[7]*av1 + Ai[8]*av2;
      h_bb[n_el_ents*ie + 3] = ab[0].min;
      h_bb[n_el_ents*ie + 4] = ab[1].min;
      h_bb[n_el_ents*ie + 5] = ab[2].min;
      h_bb[n_el_ents*ie + 6] = ab[0].max;
      h_bb[n_el_ents*ie + 7] = ab[1].max;
      h_bb[n_el_ents*ie + 8] = ab[2].max;
 
      // Expand the rotated-frame bounds by aabb_diag_len, then premultiply
      // A by the scaling matrix L = diag(di0,di1,di2) to obtain L*A. This
      // maps a physical displacement from the OBB center to [-1,1]^3.
      tb[0].min -= aabb_diag_len;
      tb[0].max += aabb_diag_len;
      tb[1].min -= aabb_diag_len;
      tb[1].max += aabb_diag_len;
      tb[2].min -= aabb_diag_len;
      tb[2].max += aabb_diag_len;
      const double di0 = 2/((1+tol)*(tb[0].max-tb[0].min)),
                   di1 = 2/((1+tol)*(tb[1].max-tb[1].min)),
                   di2 = 2/((1+tol)*(tb[2].max-tb[2].min));
      h_bb[n_el_ents*ie + 9 ]=di0*A[0];
      h_bb[n_el_ents*ie + 10]=di0*A[1];
      h_bb[n_el_ents*ie + 11]=di0*A[2];
      h_bb[n_el_ents*ie + 12]=di1*A[3];
      h_bb[n_el_ents*ie + 13]=di1*A[4];
      h_bb[n_el_ents*ie + 14]=di1*A[5];
      h_bb[n_el_ents*ie + 15]=di2*A[6];
      h_bb[n_el_ents*ie + 16]=di2*A[7];
      h_bb[n_el_ents*ie + 17]=di2*A[8];
   }
}
 
void FindPointsGSLIB::FindPointsSurfSetup3(DevStruct &devs,
                                           const double *const elx[3],
                                           const unsigned n,
                                           const uint nel,
                                           const unsigned m,
                                           const double bbox_rel_size_inc,
                                           const uint local_hash_size,
                                           const uint global_hash_size,
                                           const int rD,
                                           const Vector *aabb_sz_inc)
{
   // compute element bounding boxes.
   const bool store_obb = obb_check;
   const int sd = 3;
   const int n_box_ents = store_obb ? (3*sd + sd*sd) : (2*sd);
   devs.bb.SetSize(nel*n_box_ents);
   if (rD == 1)
   {
      obboxedge_calc_3(devs.bb, elx, n, nel, m, bbox_rel_size_inc, store_obb);
   }
   else if (rD == 2)
   {
      obboxsurf_calc_3(devs.bb, elx, n, nel, m, bbox_rel_size_inc, store_obb);
   }
   else
   {
      MFEM_ABORT("FindPointsGSLIB::FindPointsSurfSetup3: rD must be 1 or 2");
   }
 
   auto h_bb = devs.bb.HostReadWrite();
   VerifyAABBPadLayout(aabb_sz_inc, nel, sd);
   const int aabb_sz_inc_size = aabb_sz_inc ? aabb_sz_inc->Size() : 0;
 
   Vector elmin(3*nel), elmax(3*nel);
   for (uint i = 0; i < nel; i++)
   {
      const int min_off = (store_obb ? sd : 0) + n_box_ents*i;
      const int max_off = (store_obb ? 2*sd : sd) + n_box_ents*i;
      for (int d = 0; d < 3; d++)
      {
         const double pad = GetAABBPad(aabb_sz_inc, aabb_sz_inc_size,
                                       nel, sd, i, d);
         if (pad > 0.0)
         {
            h_bb[min_off + d] -= pad;
            h_bb[max_off + d] += pad;
         }
      }
      elmin(i) = h_bb[min_off + 0];
      elmin(i + nel) = h_bb[min_off + 1];
      elmin(i + 2*nel) = h_bb[min_off + 2];
      elmax(i) = h_bb[max_off + 0];
      elmax(i + nel) = h_bb[max_off + 1];
      elmax(i + 2*nel) = h_bb[max_off + 2];
   }
 
   // Skip local map construction for empty partition.
   if (nel > 0)
   {
      BBoxTensorGridMap bbmap(elmin, elmax, nel, 3, local_hash_size, true);
      devs.lh_min.HostWrite();    devs.lh_min = bbmap.GetGridMin();
      devs.lh_fac.HostWrite();    devs.lh_fac = bbmap.GetGridFac();
      devs.lh_offset.HostWrite(); devs.lh_offset = bbmap.GetGridMap();
      devs.lh_nx = bbmap.GetGridN()[0];
   }
   else
   {
      devs.lh_min.SetSize(sd);
      devs.lh_min = 0.0;
      devs.lh_fac.SetSize(sd);
      devs.lh_fac = 0.0;
      devs.lh_offset.SetSize(2);
      devs.lh_offset = 1;
      devs.lh_nx = 1;
   }
 
#ifdef MFEM_USE_MPI
   // build global map
   if (gsl_comm->np > 1)
   {
      GlobalBBoxTensorGridMap gbbmap(gsl_comm->c, elmin, elmax,
                                     nel, 3, global_hash_size, true);
      devs.gh_min = gbbmap.GetGridMin();
      devs.gh_fac = gbbmap.GetGridFac();
      devs.gh_offset = gbbmap.GetGridMap();
      devs.gh_nx = gbbmap.GetGridN()[0];
   }
#endif
}
 
void FindPointsGSLIB::FindPointsEdgeSetup2(DevStruct &devs,
                                           const double *const elx[2],
                                           const unsigned n,
                                           const uint nel,
                                           const unsigned m,
                                           const double bbox_rel_size_inc,
                                           const uint local_hash_size,
                                           const uint global_hash_size,
                                           const Vector *aabb_sz_inc)
{
   // compute element bounding boxes.
   const bool store_obb = obb_check;
   const int sd = 2;
   const int n_box_ents = store_obb ? (3*sd + sd*sd) : (2*sd);
   devs.bb.SetSize(nel*n_box_ents);
   obboxedge_calc_2(devs.bb, elx, n, nel, m, bbox_rel_size_inc, store_obb);
 
   auto h_bb = devs.bb.HostReadWrite();
   VerifyAABBPadLayout(aabb_sz_inc, nel, sd);
   const int aabb_sz_inc_size = aabb_sz_inc ? aabb_sz_inc->Size() : 0;
 
   Vector elmin(2*nel), elmax(2*nel);
   for (uint i = 0; i < nel; i++)
   {
      const int min_off = (store_obb ? sd : 0) + n_box_ents*i;
      const int max_off = (store_obb ? 2*sd : sd) + n_box_ents*i;
      for (int d = 0; d < 2 && aabb_sz_inc; d++)
      {
         const double pad = GetAABBPad(aabb_sz_inc, aabb_sz_inc_size,
                                       nel, sd, i, d);
         if (pad > 0.0)
         {
            h_bb[min_off+d] -= pad;
            h_bb[max_off+d] += pad;
         }
      }
      elmin(i) = h_bb[min_off + 0];
      elmin(i + nel) = h_bb[min_off + 1];
      elmax(i) = h_bb[max_off + 0];
      elmax(i + nel) = h_bb[max_off + 1];
   }
 
   if (nel > 0)
   {
      BBoxTensorGridMap bbmap(elmin, elmax, nel, 2, local_hash_size, true);
      devs.lh_min.HostWrite();    devs.lh_min = bbmap.GetGridMin();
      devs.lh_fac.HostWrite();    devs.lh_fac = bbmap.GetGridFac();
      devs.lh_offset.HostWrite(); devs.lh_offset = bbmap.GetGridMap();
      devs.lh_nx = bbmap.GetGridN()[0];
   }
   else
   {
      devs.lh_min.SetSize(sd);
      devs.lh_min = 0.0;
      devs.lh_fac.SetSize(sd);
      devs.lh_fac = 0.0;
      devs.lh_offset.SetSize(2);
      devs.lh_offset = 1;
      devs.lh_nx = 1;
   }
 
#ifdef MFEM_USE_MPI
   // build global map
   if (gsl_comm->np > 1)
   {
      GlobalBBoxTensorGridMap gbbmap(gsl_comm->c, elmin, elmax,
                                     nel, 2, global_hash_size, true);
      devs.gh_min = gbbmap.GetGridMin();
      devs.gh_fac = gbbmap.GetGridFac();
      devs.gh_offset = gbbmap.GetGridMap();
      devs.gh_nx = gbbmap.GetGridN()[0];
   }
#endif
}
 
// evaluate the basis function, 1st derivative, and 2nd derivative at
// x for the i-th lagrange basis function defined by the gll points z and
// lagrange coefficients
static void lagrange_eval_second_derivative(double *p0, double x, int i,
                                            const double *z,
                                            const double *lagrangeCoeff,
                                            int pN)
{
   double u0 = 1, u1 = 0, u2 = 0;
   for (int j=0; j<pN; ++j)
   {
      if (i!=j)
      {
         double d_j = 2 * (x-z[j]);
         u2 = d_j * u2 + u1;
         u1 = d_j * u1 + u0;
         u0 = d_j * u0;
      }
   }
   double *p1 = p0 + pN, *p2 = p0 + 2 * pN;
   p0[i] = lagrangeCoeff[i] * u0;
   p1[i] = 2.0 * lagrangeCoeff[i] * u1;
   p2[i] = 8.0 * lagrangeCoeff[i] * u2;
}
 
void FindPointsGSLIB::SetupSurf(Mesh &m,
                                const double bbox_rel_size_inc,
                                const double newt_tol)
{
   SetupSurfBase(m, bbox_rel_size_inc, nullptr, newt_tol);
}
 
void FindPointsGSLIB::SetupSurfWithAABBExpansion(Mesh &m,
                                                 const Vector &aabb_sz_inc,
                                                 const double newt_tol)
{
   SetupSurfBase(m, 0.0, &aabb_sz_inc, newt_tol);
}
 
void FindPointsGSLIB::SetupSurfBase(Mesh &m,
                                    const double bbox_rel_size_inc,
                                    const Vector *aabb_sz_inc,
                                    const double newt_tol)
{
   // EnsureNodes call could be useful if the mesh is 1st order and has no gridfunction defined
   MFEM_VERIFY(m.GetNodes() != nullptr, "Mesh nodes are required.");
 
   // call FreeData if FindPointsGSLIB::Setup has been called already
   if (setupflag) { FreeData(); }
   obb_check = (aabb_sz_inc == nullptr);
 
   SetupCrystal();
 
   mesh     = &m;
   dim      = mesh->Dimension();       // This is reference dimension
   spacedim = mesh->SpaceDimension();  // This is physical dimension
   MFEM_VERIFY(dim < 3, "Configuration not supported yet.");
   MFEM_VERIFY(dim < spacedim,
               "Surface setup is only for surface meshes.");
 
   bool supported_surf_elem = true;
   for (int e = 0; e < mesh->GetNE() && supported_surf_elem; e++)
   {
      const Element::Type t = mesh->GetElementType(e);
      supported_surf_elem = (t == Element::SEGMENT ||
                             t == Element::QUADRILATERAL ||
                             t == Element::TRIANGLE);
   }
   MFEM_VERIFY(supported_surf_elem,
               "FindPointsGSLIB surface mesh support: only SEGMENT, "
               "QUADRILATERAL, and TRIANGLE elements are supported.");
   MFEM_VERIFY(dim < spacedim,
               "FindPointsGSLIB surface setup only supports surface meshes.");
 
   const int meshOrder = m.GetNodes()->FESpace()->GetMaxElementOrder();
   unsigned dof1D      = meshOrder + 1;
 
   SetupSplitMeshesAndIntegrationRules(meshOrder);
   GetNodalValues(mesh->GetNodes(), gsl_mesh);
 
   mesh_points_cnt     = gsl_mesh.Size()/spacedim;
   DEV.local_hash_size = mesh_points_cnt;
   DEV.dof1d           = (int)dof1D;
   DEV.newt_tol        = newt_tol;
 
   unsigned nr = dof1D;
   unsigned mr = 2*dof1D;
   if (spacedim==2)
   {
      double * const elx[2] =
      {
         mesh_points_cnt == 0 ? nullptr : &gsl_mesh(0),
         mesh_points_cnt == 0 ? nullptr : &gsl_mesh(mesh_points_cnt)
      };
      FindPointsEdgeSetup2(DEV,
                           elx,
                           nr,
                           NE_split_total,
                           mr,
                           bbox_rel_size_inc,
                           mesh_points_cnt,
                           mesh_points_cnt,
                           aabb_sz_inc);
   }
   else if (spacedim==3)
   {
      double * const elx[3] =
      {
         mesh_points_cnt == 0 ? nullptr : &gsl_mesh(0),
         mesh_points_cnt == 0 ? nullptr : &gsl_mesh(mesh_points_cnt),
         mesh_points_cnt == 0 ? nullptr : &gsl_mesh(2*mesh_points_cnt)
      };
      FindPointsSurfSetup3(DEV,
                           elx,
                           nr,
                           NE_split_total,
                           mr,
                           bbox_rel_size_inc,
                           mesh_points_cnt,
                           mesh_points_cnt,
                           dim,
                           aabb_sz_inc);
   }
 
   // If we are applying absolute AABB expansion, compute bdr_tol such that
   // any point found within a bounding box gets marked CODE_BORDER even
   // when it is not actually on the boundary.
   if (aabb_sz_inc)
   {
      const int n_box_ents = obb_check ?
                             (3*spacedim + spacedim*spacedim) : (2*spacedim);
      const int min_off = obb_check ? spacedim : 0;
      const int max_off = obb_check ? 2*spacedim : spacedim;
      constexpr double bdr_tol_inflate = 1.01;
      double max_diag2 = bdr_tol;
      auto h_bb = DEV.bb.HostRead();
      for (int e = 0; e < NE_split_total; e++)
      {
         double diag2 = 0.0;
         for (int d = 0; d < spacedim; d++)
         {
            const double lenx = h_bb[e*n_box_ents + max_off + d] -
                                h_bb[e*n_box_ents + min_off + d];
            diag2 += lenx*lenx;
         }
         max_diag2 = std::max(max_diag2,
                              diag2 * bdr_tol_inflate * bdr_tol_inflate);
      }
#ifdef MFEM_USE_MPI
      MPI_Allreduce(MPI_IN_PLACE, &max_diag2, 1, MPI_DOUBLE, MPI_MAX,
                    gsl_comm->c);
#endif
      bdr_tol = max_diag2; // this is cube/square diagonal length squared
   }
 
   // Compute avg element size in the mesh to set surface distance tolerance.
   DEV.surf_dist_tol = 0.0;
   for (int e = 0; e < mesh->GetNE(); e++)
   {
      DEV.surf_dist_tol += mesh->GetElementVolume(e);
   }
   int nelem = NE_split_total;
#ifdef MFEM_USE_MPI
   MPI_Allreduce(MPI_IN_PLACE, &DEV.surf_dist_tol, 1, MPI_DOUBLE, MPI_SUM,
                 gsl_comm->c);
   MPI_Allreduce(MPI_IN_PLACE, &nelem, 1, MPI_INT, MPI_SUM, gsl_comm->c);
#endif
   DEV.surf_dist_tol /= nelem;
   DEV.surf_dist_tol *= 1e-16; // dist^2 tolerance so we use (1e-8)^2.
 
   // Setup gll points and lagrange coefficient data for findpoints
   Vector gll1dtemp(DEV.dof1d),
          lagcoefftemp(DEV.dof1d),
          wtendtemp(6*DEV.dof1d);
 
   gslib::lobatto_nodes(gll1dtemp.GetData(), DEV.dof1d);
   gslib::gll_lag_setup(lagcoefftemp.GetData(), DEV.dof1d);
 
   for (int i=0; i<DEV.dof1d; i++)
   {
      lagrange_eval_second_derivative(wtendtemp.GetData(), -1.0, i,
                                      gll1dtemp.GetData(),
                                      lagcoefftemp.GetData(), DEV.dof1d);
      lagrange_eval_second_derivative(wtendtemp.GetData()+3*DEV.dof1d, 1.0, i,
                                      gll1dtemp.GetData(),
                                      lagcoefftemp.GetData(), DEV.dof1d);
   }
 
   DEV.wtend.SetSize(6*DEV.dof1d);
   auto h_wtend = DEV.wtend.HostWrite();
   for (int i = 0; i < 6*DEV.dof1d; i++) { h_wtend[i] = wtendtemp[i]; }
 
   DEV.gll1d.SetSize(DEV.dof1d);
   auto h_gll1d = DEV.gll1d.HostWrite();
   for (int i = 0; i < DEV.dof1d; i++) { h_gll1d[i] = gll1dtemp[i]; }
 
   DEV.lagcoeff.SetSize(DEV.dof1d);
   auto h_lagcoeff = DEV.lagcoeff.HostWrite();
   for (int i = 0; i < DEV.dof1d; i++) { h_lagcoeff[i] = lagcoefftemp[i]; }
 
   DEV.cr = cr;
   setupflag = true;
}
 
void FindPointsGSLIB::FindPoints(const Vector &point_pos,
                                 const int point_pos_ordering)
{
   MFEM_VERIFY(setupflag, "Use FindPointsGSLIB::Setup before finding points.");
   if (dim != spacedim)
   {
      FindPointsSurf(point_pos, point_pos_ordering);
      return;
   }
   bool dev_mode = (point_pos.UseDevice() && Device::IsEnabled());
   points_cnt = point_pos.Size() / dim;
   gsl_code.SetSize(points_cnt);
   gsl_proc.SetSize(points_cnt);
   gsl_elem.SetSize(points_cnt);
   gsl_ref.SetSize(points_cnt * dim);
   gsl_dist.SetSize(points_cnt);
 
   bool tensor_product_only = mesh->GetNE() == 0 ||
                              (mesh->GetNumGeometries(dim) == 1 &&
                               (mesh->GetElementType(0)==Element::QUADRILATERAL ||
                                mesh->GetElementType(0) == Element::HEXAHEDRON));
#ifdef MFEM_USE_MPI
   MPI_Allreduce(MPI_IN_PLACE, &tensor_product_only, 1, MFEM_MPI_CXX_BOOL,
                 MPI_LAND, gsl_comm->c);
#endif
 
   if (dev_mode && tensor_product_only)
   {
#if GSLIB_RELEASE_VERSION == 10007
      if (!gpu_to_cpu_fallback)
      {
         MFEM_ABORT("Either update to gslib v1.0.9 for GPU support "
                    "or use SetGPUtoCPUFallback to use host-functions. See "
                    "INSTALL for instructions to update GSLIB.");
      }
#else
      FindPointsOnDevice(point_pos, point_pos_ordering);
      return;
#endif
   }
 
   auto pp = point_pos.HostRead();
   auto xvFill = [&](const double *xv_base[], unsigned xv_stride[])
   {
      for (int d = 0; d < dim; d++)
      {
         if (point_pos_ordering == Ordering::byNODES)
         {
            xv_base[d] = pp + d*points_cnt;
            xv_stride[d] = sizeof(double);
         }
         else
         {
            xv_base[d] = pp + d;
            xv_stride[d] = dim*sizeof(double);
         }
      }
   };
 
   if (dim == 2)
   {
      auto *findptsData = (gslib::findpts_data_2 *)this->fdataD;
      const double *xv_base[2];
      unsigned xv_stride[2];
      xvFill(xv_base, xv_stride);
      findpts_2(gsl_code.GetData(), sizeof(unsigned int),
                gsl_proc.GetData(), sizeof(unsigned int),
                gsl_elem.GetData(), sizeof(unsigned int),
                gsl_ref.GetData(),  sizeof(double) * dim,
                gsl_dist.GetData(), sizeof(double),
                xv_base, xv_stride, points_cnt, findptsData);
   }
   else  // dim == 3
   {
      auto *findptsData = (gslib::findpts_data_3 *)this->fdataD;
      const double *xv_base[3];
      unsigned xv_stride[3];
      xvFill(xv_base, xv_stride);
      findpts_3(gsl_code.GetData(), sizeof(unsigned int),
                gsl_proc.GetData(), sizeof(unsigned int),
                gsl_elem.GetData(), sizeof(unsigned int),
                gsl_ref.GetData(),  sizeof(double) * dim,
                gsl_dist.GetData(), sizeof(double),
                xv_base, xv_stride, points_cnt,
                findptsData);
   }
 
   // Set the element number and reference position to 0 for points not found
   for (int i = 0; i < points_cnt; i++)
   {
      if (gsl_code[i] == 2 ||
          (gsl_code[i] == 1 && gsl_dist(i) > bdr_tol))
      {
         gsl_elem[i] = 0;
         for (int d = 0; d < dim; d++) { gsl_ref(i*dim + d) = -1.; }
         gsl_code[i] = 2;
         gsl_proc[i] = gsl_comm->id;
      }
   }
 
   // Map element number for simplices, and ref_pos from [-1,1] to [0,1] for
   // both simplices and quads. Also sets code to 1 for points found on element
   // faces/edges.
   MapRefPosAndElemIndices();
}
 
#if GSLIB_RELEASE_VERSION >= 10009
 
slong lfloor(double x) { return floor(x); }
 
// Local hash mesh index in 1D for a given point
ulong hash_index_1(double low, double fac, ulong n, double x)
{
   const slong i = lfloor((x - low) * fac);
   return i < 0 ? 0 : (n - 1 < (ulong)i ? n - 1 : (ulong)i);
}
 
// Local hash mesh index in 3D for a given point
ulong hash_index_3(const gslib::hash_data_3 *p, const double x[3])
{
   const ulong n = p->hash_n;
   return (hash_index_1(p->bnd[2].min, p->fac[2], n, x[2]) * n +
           hash_index_1(p->bnd[1].min, p->fac[1], n, x[1])) *
          n +
          hash_index_1(p->bnd[0].min, p->fac[0], n, x[0]);
}
 
// Local hash mesh index in 2D for a given point
ulong hash_index_2(const gslib::hash_data_2 *p, const double x[2])
{
   const ulong n = p->hash_n;
   return (hash_index_1(p->bnd[1].min, p->fac[1], n, x[1])) * n
          + hash_index_1(p->bnd[0].min, p->fac[0], n, x[0]);
}
 
void FindPointsGSLIB::SetupDevice()
{
   auto *findptsData3 = (gslib::findpts_data_3 *)this->fdataD;
   auto *findptsData2 = (gslib::findpts_data_2 *)this->fdataD;
 
   DEV.newt_tol = dim == 2 ? findptsData2->local.tol : findptsData3->local.tol;
   if (dim == 3)
   {
      DEV.hash3 = &findptsData3->hash;
   }
   else
   {
      DEV.hash2 = &findptsData2->hash;
   }
   DEV.cr = dim == 2 ? &findptsData2->cr : &findptsData3->cr;
 
   gsl_mesh.UseDevice(true);
 
   int n_box_ents = 3*dim + dim*dim;
   DEV.bb.UseDevice(true); DEV.bb.SetSize(n_box_ents*NE_split_total);
   auto p_bb = DEV.bb.HostWrite();
   Vector elmin(dim*NE_split_total), elmax(dim*NE_split_total);
 
   const int dim2 = dim*dim;
   if (dim == 3)
   {
      for (int e = 0; e < NE_split_total; e++)
      {
         auto box = findptsData3->local.obb[e];
         for (int d = 0; d < dim; d++)
         {
            p_bb[n_box_ents*e + d] = box.c0[d];
            p_bb[n_box_ents*e + dim + d] = box.x[d].min;
            p_bb[n_box_ents*e + 2*dim + d] = box.x[d].max;
            elmin(d*NE_split_total + e) = box.x[d].min;
            elmax(d*NE_split_total + e) = box.x[d].max;
         }
         for (int d = 0; d < dim2; ++d)
         {
            p_bb[n_box_ents*e + 3*dim + d] = box.A[d];
         }
      }
   }
   else
   {
      for (int e = 0; e < NE_split_total; e++)
      {
         auto box = findptsData2->local.obb[e];
         for (int d = 0; d < dim; d++)
         {
            p_bb[n_box_ents*e + d] = box.c0[d];
            p_bb[n_box_ents*e + dim + d] = box.x[d].min;
            p_bb[n_box_ents*e + 2*dim + d] = box.x[d].max;
            elmin(d*NE_split_total + e) = box.x[d].min;
            elmax(d*NE_split_total + e) = box.x[d].max;
         }
         for (int d = 0; d < dim2; ++d)
         {
            p_bb[n_box_ents*e + 3*dim + d] = box.A[d];
         }
      }
   }
 
   DEV.lh_min.UseDevice(true); DEV.lh_min.SetSize(dim);
   DEV.lh_fac.UseDevice(true); DEV.lh_fac.SetSize(dim);
   if (dim == 2)
   {
      auto hash = findptsData2->local.hd;
      auto p_loc_hash_min = DEV.lh_min.HostWrite();
      auto p_loc_hash_fac = DEV.lh_fac.HostWrite();
      for (int d = 0; d < dim; d++)
      {
         p_loc_hash_min[d] = hash.bnd[d].min;
         p_loc_hash_fac[d] = hash.fac[d];
      }
      DEV.lh_nx = hash.hash_n;
   }
   else
   {
      auto hash = findptsData3->local.hd;
      auto p_loc_hash_min = DEV.lh_min.HostWrite();
      auto p_loc_hash_fac = DEV.lh_fac.HostWrite();
      for (int d = 0; d < dim; d++)
      {
         p_loc_hash_min[d] = hash.bnd[d].min;
         p_loc_hash_fac[d] = hash.fac[d];
      }
      DEV.lh_nx = hash.hash_n;
   }
 
   int h_o_size = dim == 2 ?
                  findptsData2->local.hd.offset[(int)std::pow(DEV.lh_nx, dim)] :
                  findptsData3->local.hd.offset[(int)std::pow(DEV.lh_nx, dim)];
 
   DEV.lh_offset.SetSize(h_o_size);
   auto p_ou_offset = DEV.lh_offset.HostWrite();
   for (int i = 0; i < h_o_size; i++)
   {
      p_ou_offset[i] = dim == 2 ? findptsData2->local.hd.offset[i] :
                       findptsData3->local.hd.offset[i];
   }
 
   DEV.wtend.UseDevice(true);
   DEV.wtend.SetSize(6*DEV.dof1d);
   auto *h_wtend = DEV.wtend.HostWrite();
   double *src_wtend = dim == 2 ? findptsData2->local.fed.wtend[0] :
                       findptsData3->local.fed.wtend[0];
   for (int i = 0; i < 6*DEV.dof1d; i++) { h_wtend[i] = src_wtend[i]; }
 
   // Get gll points
   DEV.gll1d.UseDevice(true);
   DEV.gll1d.SetSize(DEV.dof1d);
   auto *h_gll1d = DEV.gll1d.HostWrite();
   double *src_gll1d = dim == 2 ? findptsData2->local.fed.z[0] :
                       findptsData3->local.fed.z[0];
   for (int i = 0; i < DEV.dof1d; i++) { h_gll1d[i] = src_gll1d[i]; }
 
   DEV.lagcoeff.UseDevice(true);
   DEV.lagcoeff.SetSize(DEV.dof1d);
   auto *h_lagcoeff = DEV.lagcoeff.HostWrite();
   double *src_lagcoeff = dim == 2 ? findptsData2->local.fed.lag_data[0] :
                          findptsData3->local.fed.lag_data[0];
   for (int i = 0; i < DEV.dof1d; i++) { h_lagcoeff[i] = src_lagcoeff[i]; }
 
   DEV.setup_device = true;
}
 
void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
                                         const int point_pos_ordering)
{
   if (!DEV.setup_device)
   {
      SetupDevice();
   }
   DEV.find_device = true;
 
   const int id = gsl_comm->id, np = gsl_comm->np;
   const int ddim = dim;
 
   gsl_mfem_ref.SetSize(points_cnt * dim);
   gsl_mfem_elem.SetSize(points_cnt);
   gsl_ref.UseDevice(true);
   gsl_dist.UseDevice(true);
 
   // Initialize proc to id for all points. Code and dist will be done
   // inside kernels.
   gsl_proc.HostWrite();
   gsl_proc = id;
 
   if (dim == 2)
   {
      FindPointsLocal2(point_pos, point_pos_ordering, gsl_code, gsl_elem,
                       gsl_ref, gsl_dist, points_cnt);
   }
   else
   {
      FindPointsLocal3(point_pos, point_pos_ordering, gsl_code, gsl_elem,
                       gsl_ref, gsl_dist, points_cnt);
   }
 
   // Tolerance for point to be marked as on element edge/face based on the
   // obtained reference-space coordinates.
   double rbtol = 1e-12; // must match MapRefPosAndElemIndices for consistency
 
   if (np == 1)
   {
      auto d_gsl_code = gsl_code.ReadWrite(); // no memory transfer
      auto d_gsl_ref = gsl_ref.Read(); // no memory transfer
      auto d_gsl_dist = gsl_dist.Read(); // no memory transfer
      auto d_gsl_elem = gsl_elem.Read(); // no memory transfer
      auto d_gsl_mfem_ref = gsl_mfem_ref.Write();
      auto d_gsl_mfem_elem = gsl_mfem_elem.Write();
 
      const int pts_cnt = points_cnt;
      const double bdr_t = bdr_tol;
 
      FindPointsDeviceSetCode(pts_cnt, ddim, rbtol, bdr_t, d_gsl_code,
                              d_gsl_ref, d_gsl_dist, d_gsl_elem,
                              d_gsl_mfem_ref, d_gsl_mfem_elem);
      return;
   }
 
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
   /* send unfound and border points to global hash cells */
   struct gslib::array hash_pt, src_pt, out_pt;
 
   struct srcPt_t
   {
      double x[3];
      unsigned int index, proc;
   };
 
   struct outPt_t
   {
      double r[3], dist2;
      unsigned int index, code, el, proc;
   };
 
   {
      // Sync from device to host to send points to other ranks
      auto h_gsl_code = gsl_code.HostReadWrite(); // read now, write later
      auto h_pp = point_pos.HostRead();
 
      int index;
      struct srcPt_t *pt;
 
      array_init(struct srcPt_t, &hash_pt, points_cnt);
      pt = (struct srcPt_t *)hash_pt.ptr;
 
      auto x = new double[dim];
      for (index = 0; index < points_cnt; ++index)
      {
         const int code_i = h_gsl_code[index];
         if (code_i != CODE_INTERNAL)
         {
            for (int d = 0; d < dim; ++d)
            {
               int idx = point_pos_ordering == 0 ?
                         index + d*points_cnt :
                         index*dim + d;
               x[d] = h_pp[idx];
            }
            const auto hi = dim == 2 ? hash_index_2(DEV.hash2, x) :
                            hash_index_3(DEV.hash3, x);
            for (int d = 0; d < dim; ++d)
            {
               pt->x[d] = x[d];
            }
            pt->index = index;
            pt->proc = hi % np;
            ++pt;
         }
      }
      delete[] x;
      hash_pt.n = pt - (struct srcPt_t *)hash_pt.ptr;
      sarray_transfer(struct srcPt_t, &hash_pt, proc, 1, DEV.cr);
   }
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
 
   /* look up points in hash cells, route to possible procs */
   {
      const unsigned int *const hash_offset = dim == 2 ? DEV.hash2->offset :
                                              DEV.hash3->offset;
      int count = 0;
      unsigned int *proc, *proc_p;
      const struct srcPt_t *p = (struct srcPt_t *)hash_pt.ptr,
                            *const pe = p + hash_pt.n;
      struct srcPt_t *q;
 
      for (; p != pe; ++p)
      {
         const int hi = dim == 2 ? hash_index_2(DEV.hash2, p->x)/np :
                        hash_index_3(DEV.hash3, p->x)/np;
         const int i = hash_offset[hi], ie = hash_offset[hi + 1];
         count += ie - i;
      }
 
      Array<unsigned int> proc_array(count);
      proc = proc_array.GetData();
      proc_p = proc;
      array_init(struct srcPt_t, &src_pt, count);
      q = (struct srcPt_t *)src_pt.ptr;
 
      p = (struct srcPt_t *)hash_pt.ptr;
      for (; p != pe; ++p)
      {
         const int hi = dim == 2 ? hash_index_2(DEV.hash2, p->x)/np :
                        hash_index_3(DEV.hash3, p->x)/np;
         int i = hash_offset[hi];
         const int ie = hash_offset[hi + 1];
         for (; i != ie; ++i)
         {
            const unsigned int pp = hash_offset[i];
            /* don't send back to where it just came from */
            if (pp == p->proc)
            {
               continue;
            }
            *proc_p++ = pp;
            *q++ = *p;
         }
      }
 
      array_free(&hash_pt);
      src_pt.n = proc_p - proc;
 
      sarray_transfer_ext(struct srcPt_t, &src_pt, proc, sizeof(uint), DEV.cr);
   }
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
 
   /* look for other procs' points, send back */
   {
      int n = src_pt.n;
      const struct srcPt_t *spt;
      struct outPt_t *opt;
      array_init(struct outPt_t, &out_pt, n);
      out_pt.n = n;
      spt = (struct srcPt_t *)src_pt.ptr;
      opt = (struct outPt_t *)out_pt.ptr;
      for (; n; --n, ++spt, ++opt)
      {
         opt->index = spt->index;
         opt->proc = spt->proc;
      }
      spt = (struct srcPt_t *)src_pt.ptr;
      opt = (struct outPt_t *)out_pt.ptr;
 
      n = out_pt.n;
      Vector gsl_ref_l, gsl_dist_l;
      gsl_ref_l.UseDevice(true); gsl_ref_l.SetSize(n*dim);
      gsl_dist_l.UseDevice(true); gsl_dist_l.SetSize(n);
 
      Vector point_pos_l;
      point_pos_l.UseDevice(true); point_pos_l.SetSize(n*dim);
      auto h_ptr_ppl = point_pos_l.HostWrite();
 
      Array<unsigned int> gsl_code_l(n), gsl_elem_l(n);
 
      for (int point = 0; point < n; ++point)
      {
         for (int d = 0; d < dim; d++)
         {
            int idx = point_pos_ordering == 0 ? point + d*n : point*dim + d;
            h_ptr_ppl[idx] = spt[point].x[d];
         }
      }
 
      if (dim == 2)
      {
         FindPointsLocal2(point_pos_l, point_pos_ordering, gsl_code_l,
                          gsl_elem_l, gsl_ref_l, gsl_dist_l, n);
      }
      else
      {
         FindPointsLocal3(point_pos_l, point_pos_ordering, gsl_code_l,
                          gsl_elem_l, gsl_ref_l, gsl_dist_l, n);
      }
 
      gsl_ref_l.HostRead();
      gsl_dist_l.HostRead();
      gsl_code_l.HostRead();
      gsl_elem_l.HostRead();
 
      // unpack arrays into opt
      for (int point = 0; point < n; point++)
      {
         opt[point].code = AsConst(gsl_code_l)[point];
         if (opt[point].code == CODE_NOT_FOUND)
         {
            continue;
         }
         opt[point].el   = AsConst(gsl_elem_l)[point];
         opt[point].dist2 = AsConst(gsl_dist_l)[point];
         for (int d = 0; d < dim; ++d)
         {
            opt[point].r[d] = AsConst(gsl_ref_l)[dim * point + d];
         }
         // for found points set gsl_code using reference space coords.
         IntegrationPoint ip;
         if (dim == 2)
         {
            ip.Set2(0.5*opt[point].r[0]+0.5, 0.5*opt[point].r[1]+0.5);
         }
         else
         {
            ip.Set3(0.5*opt[point].r[0]+0.5, 0.5*opt[point].r[1]+0.5,
                    0.5*opt[point].r[2]+0.5);
         }
         const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(opt[point].el);
         const Geometry::Type gt = fe->GetGeomType();
         int setcode = Geometry::CheckPoint(gt, ip, -rbtol) ?
                       CODE_INTERNAL : CODE_BORDER;
         opt[point].code = setcode==CODE_BORDER && opt[point].dist2>bdr_tol ?
                           CODE_NOT_FOUND : setcode;
      }
 
      array_free(&src_pt);
 
      /* group by code to eliminate unfound points */
      sarray_sort(struct outPt_t, opt, out_pt.n, code, 0, &DEV.cr->data);
 
      n = out_pt.n;
      while (n && opt[n - 1].code == CODE_NOT_FOUND)
      {
         --n;
      }
      out_pt.n = n;
 
      sarray_transfer(struct outPt_t, &out_pt, proc, 1, DEV.cr);
   }
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
 
   auto h_gsl_code = gsl_code.HostReadWrite();
   auto h_gsl_dist = gsl_dist.HostReadWrite();
   auto h_gsl_ref = gsl_ref.HostReadWrite();
   auto h_gsl_mfem_ref = gsl_mfem_ref.HostReadWrite();
   auto h_gsl_elem = gsl_elem.HostReadWrite();
   auto h_gsl_mfem_elem = gsl_mfem_elem.HostReadWrite();
   auto h_gsl_proc = gsl_proc.HostReadWrite();
 
   /* merge remote results with user data */
   {
      int n = out_pt.n;
      struct outPt_t *opt = (struct outPt_t *)out_pt.ptr;
      for (int i = 0; i < n; i++)
      {
         const int index = opt[i].index;
         if (h_gsl_code[index] == CODE_INTERNAL)
         {
            continue;
         }
         if (h_gsl_code[index] == CODE_NOT_FOUND ||
             opt[i].code == CODE_INTERNAL ||
             opt[i].dist2 < h_gsl_dist[index])
         {
            for (int d = 0; d < dim; ++d)
            {
               real_t rv = opt[i].r[d];
               h_gsl_ref[dim * index + d] = rv;
               h_gsl_mfem_ref[dim*index + d] = 0.5*(rv + 1.);
            }
            h_gsl_dist[index] = opt[i].dist2;
            h_gsl_proc[index] = opt[i].proc;
            h_gsl_elem[index] = opt[i].el;
            h_gsl_mfem_elem[index]   = opt[i].el;
            h_gsl_code[index] = opt[i].code;
         }
      }
      array_free(&out_pt);
   }
 
   // For points found locally, we set gsl_mfem_elem, gsl_mfem_ref, and
   // gsl_code since it was not set until now.
   for (int index = 0; index < points_cnt; index++)
   {
      if (h_gsl_code[index] == CODE_NOT_FOUND || h_gsl_proc[index] != id)
      {
         continue;
      }
      h_gsl_mfem_elem[index] = h_gsl_elem[index];
 
      bool internal = true;
      for (int k = 0; k < dim; k++)
      {
         double r_val = h_gsl_ref[index * dim + k];
         double val = 0.5 * (r_val + 1.0);
         h_gsl_mfem_ref[index * dim + k] = val;
         if (val < rbtol || val > 1.0 - rbtol)
         {
            internal = false;
         }
      }
 
      int setcode = internal ? CODE_INTERNAL : CODE_BORDER;
      h_gsl_code[index] = (setcode == CODE_BORDER &&
                           h_gsl_dist[index] > bdr_tol)
                          ? CODE_NOT_FOUND
                          : setcode;
   }
}
 
struct evalSrcPt_t
{
   double r[3];
   unsigned int index, proc, el;
};
 
struct evalOutPt_t
{
   double out;
   unsigned int index, proc;
};
 
void FindPointsGSLIB::InterpolateOnDevice(const Vector &field_in_evec,
                                          Vector &field_out,
                                          const int nel,
                                          const int ncomp,
                                          const int dof1Dsol,
                                          const int ordering)
{
   field_out.UseDevice(true);
   field_out.SetSize(points_cnt*ncomp);
   field_out = default_interp_value;
 
   DEV.dof1d_sol =  dof1Dsol;
   DEV.gll1d_sol.UseDevice(true);  DEV.gll1d_sol.SetSize(dof1Dsol);
   DEV.lagcoeff_sol.UseDevice(true); DEV.lagcoeff_sol.SetSize(dof1Dsol);
   if (DEV.dof1d_sol != DEV.dof1d || !DEV.find_device)
   {
      gslib::lobatto_nodes(DEV.gll1d_sol.HostWrite(), dof1Dsol);
      gslib::gll_lag_setup(DEV.lagcoeff_sol.HostWrite(), dof1Dsol);
   }
   else
   {
      DEV.gll1d_sol = DEV.gll1d.HostRead();
      DEV.lagcoeff_sol = DEV.lagcoeff.HostRead();
   }
 
   // field_out.HostReadWrite(); //Reads in default value from device
 
   struct gslib::array src, outpt;
   int nlocal = 0;
   /* weed out unfound points, send out */
   Array<int> gsl_elem_temp;
   Vector gsl_ref_temp;
   Array<int> index_temp;
   {
      int index;
      const unsigned int *h_code = gsl_code.HostRead(),
                          *h_proc = gsl_proc.HostRead(),
                           *h_el   = gsl_elem.HostRead();
      const double *h_r = gsl_ref.HostRead();
 
      int numSend = 0;
 
      for (index = 0; index < points_cnt; ++index)
      {
         numSend += (h_code[index] != CODE_NOT_FOUND &&
                     h_proc[index] != gsl_comm->id);
         nlocal += (h_code[index] != CODE_NOT_FOUND &&
                    h_proc[index] == gsl_comm->id);
      }
 
      gsl_elem_temp.SetSize(nlocal);
      gsl_elem_temp.HostWrite();
 
      gsl_ref_temp.SetSize(nlocal*dim);
      gsl_ref_temp.UseDevice(true);
      gsl_ref_temp.HostWrite();
 
      index_temp.SetSize(nlocal);
 
      evalSrcPt_t *pt;
      array_init(evalSrcPt_t, &src, numSend);
      pt = (evalSrcPt_t *)src.ptr;
 
      int ctr = 0;
      for (index = 0; index < points_cnt; ++index)
      {
         if (h_code[index] != CODE_NOT_FOUND && h_proc[index] != gsl_comm->id)
         {
            for (int d = 0; d < dim; ++d)
            {
               pt->r[d] = h_r[index*dim+d];
            }
            pt->index = index;
            pt->proc = h_proc[index];
            pt->el = h_el[index];
            ++pt;
         }
         else if (h_code[index] != CODE_NOT_FOUND &&
                  h_proc[index] == gsl_comm->id)
         {
            gsl_elem_temp[ctr] = h_el[index];
            for (int d = 0; d < dim; ++d)
            {
               gsl_ref_temp(dim*ctr+d) = h_r[index*dim+d];
            }
            index_temp[ctr] = index;
            ctr++;
         }
      }
 
      src.n = pt - (evalSrcPt_t *)src.ptr;
      sarray_transfer(evalSrcPt_t, &src, proc, 1, cr);
   }
 
   //evaluate points that are already local
   {
      Vector interp_vals(nlocal*ncomp);
      interp_vals.UseDevice(true);
 
      if (dim == 2)
      {
         InterpolateLocal2(field_in_evec,
                           gsl_elem_temp,
                           gsl_ref_temp,
                           interp_vals,
                           nlocal, ncomp,
                           dof1Dsol);
      }
      else
      {
         InterpolateLocal3(field_in_evec,
                           gsl_elem_temp,
                           gsl_ref_temp,
                           interp_vals,
                           nlocal, ncomp,
                           dof1Dsol);
 
      }
#ifdef MFEM_USE_MPI
      MPI_Barrier(gsl_comm->c);
#endif
 
      auto d_interp_vals = interp_vals.Read();
      auto d_index_temp  = index_temp.Read();
      auto d_field_out   = field_out.ReadWrite(); // no-op, already on device.
 
      const int interp_Offset = interp_vals.Size()/ncomp;
      const int pts_cnt = points_cnt;
 
      InterpolateDeviceScatter(nlocal, d_index_temp, d_interp_vals,
                               d_field_out, interp_Offset, ncomp, pts_cnt,
                               ordering);
   }
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
 
   if (gsl_comm->np == 1)
   {
      array_free(&src);
      return;
   }
 
   // evaluate points received from other ranks
   {
      int n = src.n;
      const evalSrcPt_t *spt;
      spt = (evalSrcPt_t *)src.ptr;
 
      // Copy to host vector
      gsl_elem_temp.SetSize(n);
      gsl_elem_temp.HostWrite();
 
      gsl_ref_temp.SetSize(n*dim);
      gsl_ref_temp.HostWrite();
 
      spt = (evalSrcPt_t *)src.ptr;
      for (int i = 0; i < n; i++, ++spt)
      {
         gsl_elem_temp[i] = spt->el;
         for (int d = 0; d < dim; d++)
         {
            gsl_ref_temp(i*dim + d) = spt->r[d];
         }
      }
 
      Vector interp_vals(n*ncomp);
      interp_vals.UseDevice(true);
      if (dim == 2)
      {
         InterpolateLocal2(field_in_evec,
                           gsl_elem_temp,
                           gsl_ref_temp,
                           interp_vals, n, ncomp, dof1Dsol);
      }
      else
      {
         InterpolateLocal3(field_in_evec,
                           gsl_elem_temp,
                           gsl_ref_temp,
                           interp_vals, n, ncomp, dof1Dsol);
      }
#ifdef MFEM_USE_MPI
      MPI_Barrier(gsl_comm->c);
#endif
      auto h_interp_vals = interp_vals.HostRead();
      auto h_field_out   = field_out.HostReadWrite();
 
      // Now the interpolated values need to be sent back component wise
      int Offset = interp_vals.Size()/ncomp;
      for (int i = 0; i < ncomp; i++)
      {
         spt = (evalSrcPt_t *)src.ptr;
         array_init(evalOutPt_t, &outpt, n);
         outpt.n = n;
         evalOutPt_t *opt = (evalOutPt_t *)outpt.ptr;
 
         for (int j = 0; j < n; j++)
         {
            opt->index = spt->index;
            opt->proc = spt->proc;
            opt->out = h_interp_vals[j + Offset*i];
            spt++;
            opt++;
         }
 
         sarray_transfer(struct evalOutPt_t, &outpt, proc, 1, cr);
 
         opt = (evalOutPt_t *)outpt.ptr;
         for (size_t index = 0; index < outpt.n; index++)
         {
            int idx = ordering == Ordering::byNODES ?
                      opt->index + i*points_cnt :
                      opt->index*ncomp + i;
            h_field_out[idx] = opt->out;
            ++opt;
         }
         array_free(&outpt);
      }
      array_free(&src);
   }
   field_out.ReadWrite();
   //finished evaluating points received from other processors.
}
 
// Hash mesh index in 3D for a given point
ulong hash_index_nd(const Vector &hash_min, const Vector &hash_fac,
                    const int n,  const Vector &x)
{
   const int dim = x.Size();
   if (dim == 2)
   {
      return (hash_index_1(hash_min(1), hash_fac(1), n, x[1])) * n
             + hash_index_1(hash_min(0), hash_fac(0), n, x[0]);
   }
   else if (dim == 3)
   {
      return (hash_index_1(hash_min(2), hash_fac(2), n, x[2]) * n +
              hash_index_1(hash_min(1), hash_fac(1), n, x[1])) * n +
             hash_index_1(hash_min(0), hash_fac(0), n, x[0]);
   }
   else
   {
      MFEM_ABORT("hash_index_nd only supports 2D and 3D cases.");
   }
}
 
void FindPointsGSLIB::FindPointsSurf(const Vector &point_pos,
                                     int point_pos_ordering)
{
   MFEM_VERIFY(setupflag, "Use FindPointsGSLIB::Setup before finding points.");
   MFEM_VERIFY(dim < spacedim, "FindPointsSurf is only for surface meshes.");
 
   bool supported_surf_elem = true; // future-proof
   for (int e = 0; e < mesh->GetNE(); e++)
   {
      const Element::Type t = mesh->GetElementType(e);
      if (t != Element::TRIANGLE &&
          t != Element::SEGMENT &&
          t != Element::QUADRILATERAL)
      {
         supported_surf_elem = false;
      }
   }
#ifdef MFEM_USE_MPI
   MPI_Allreduce(MPI_IN_PLACE, &supported_surf_elem, 1, MFEM_MPI_CXX_BOOL,
                 MPI_LAND, gsl_comm->c);
#endif
   MFEM_VERIFY(supported_surf_elem,
               "FindPointsGSLIB surface mesh support: only SEGMENT, "
               "QUADRILATERAL, and TRIANGLE elements are supported.");
 
   points_cnt = point_pos.Size()/spacedim;
 
   gsl_code.SetSize(points_cnt);
   gsl_proc.SetSize(points_cnt);
   gsl_elem.SetSize(points_cnt);
   gsl_ref.SetSize(points_cnt*dim);
   gsl_dist.SetSize(points_cnt);
   gsl_mfem_ref.SetSize(points_cnt * dim);
   gsl_mfem_elem.SetSize(points_cnt);
 
   bool use_dev = point_pos.UseDevice();
 
   gsl_ref.UseDevice(use_dev);
   gsl_dist.UseDevice(use_dev);
 
   if (spacedim==2)
   {
      FindPointsEdgeLocal2(point_pos,
                           point_pos_ordering,
                           gsl_code,
                           gsl_elem,
                           gsl_ref,
                           gsl_dist,
                           points_cnt);
   }
   else
   {
      if (dim == 1)
      {
         FindPointsEdgeLocal3(point_pos,
                              point_pos_ordering,
                              gsl_code,
                              gsl_elem,
                              gsl_ref,
                              gsl_dist,
                              points_cnt);
      }
      else if (dim == 2)
      {
         FindPointsSurfLocal3(point_pos,
                              point_pos_ordering,
                              gsl_code,
                              gsl_elem,
                              gsl_ref,
                              gsl_dist,
                              points_cnt);
 
      }
   }
 
   // tolerance for point to be marked as on element edge/face
   const unsigned int id = gsl_comm->id,
                      np = gsl_comm->np;
   gsl_proc.HostWrite();
   for (int i=0; i<points_cnt; i++)
   {
      gsl_proc[i] = id;
   }
   // Tolerance for point to be marked as on element edge/face based on the
   // obtained reference-space coordinates.
   const double rbtol = 1e-12; // must match MapRefPosAndElemIndices
 
   auto d_gsl_code     = gsl_code.ReadWrite(use_dev);
   auto d_gsl_ref      = gsl_ref.Read(use_dev);
   auto d_gsl_elem     = gsl_elem.Read(use_dev);
   auto d_gsl_mfem_ref = gsl_mfem_ref.Write(use_dev);
   auto d_gsl_mfem_elem = gsl_mfem_elem.Write(use_dev);
   auto d_gsl_dist     = gsl_dist.Read(use_dev);
   const int ddim = dim;
   const double dbdr_tol = bdr_tol;
   auto d_split_elem_map  = split_element_map.Read(use_dev);
   auto d_split_elem_idx  = split_element_index.Read(use_dev);
   auto d_split_elem_geom = split_element_geom.Read(use_dev);
 
   // Note: we check if the point is on element border and mark it as such.
   // We do not mark points as CODE_INTERNAL because the found solution could
   // be interior to the element even when the point is not on the surface.
   // This case is handled in the kernels.
   MFEM_FORALL(index, points_cnt,
   {
      if (d_gsl_code[index] == CODE_NOT_FOUND) { return; }
 
      if (ddim == 1)
      {
         // Segment surface: no splitting, direct copy.
         d_gsl_mfem_elem[index] = d_gsl_elem[index];
         const real_t ipx = 0.5*(d_gsl_ref[index] + 1.0);
         d_gsl_mfem_ref[index] = ipx;
         if (ipx < rbtol || ipx > 1.0 - rbtol)
         {
            d_gsl_code[index] = CODE_BORDER;
         }
      }
      else if (ddim == 2)
      {
         const int loc_id = d_gsl_elem[index];
         d_gsl_mfem_elem[index] = d_split_elem_map[loc_id];
 
         if (d_split_elem_geom[loc_id] == (int)Geometry::TRIANGLE)
         {
            // Each original triangle is split into 3 quads.
            const int tri_id = d_split_elem_idx[loc_id]; // 0, 1, or 2
 
            const double u = 0.5*(d_gsl_ref[index*2 + 0] + 1.0);
            const double v = 0.5*(d_gsl_ref[index*2 + 1] + 1.0);
            double tx, ty;
            MapSplitTriangleQuadToTriangle(tri_id, u, v, tx, ty);
 
            d_gsl_mfem_ref[index*2 + 0] = tx;
            d_gsl_mfem_ref[index*2 + 1] = ty;
 
            // Triangle edges
            if (tx < rbtol || ty < rbtol || tx + ty > 1.0 - rbtol)
            {
               d_gsl_code[index] = CODE_BORDER;
            }
         }
         else
         {
            // Quad/segment boundary/end point
            const real_t ipx = 0.5*(d_gsl_ref[index*2 + 0] + 1.0);
            const real_t ipy = 0.5*(d_gsl_ref[index*2 + 1] + 1.0);
            d_gsl_mfem_ref[index*2 + 0] = ipx;
            d_gsl_mfem_ref[index*2 + 1] = ipy;
            if (ipx < rbtol || ipx > 1.0 - rbtol ||
                ipy < rbtol || ipy > 1.0 - rbtol)
            {
               d_gsl_code[index] = CODE_BORDER;
            }
         }
      }
      if (d_gsl_code[index] == CODE_BORDER && d_gsl_dist[index] > dbdr_tol)
      {
         d_gsl_code[index] = CODE_NOT_FOUND;
      }
   });
   if (np == 1) { return; }
 
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
 
   /* send unfound and border points to global hash cells */
   struct gslib::array hash_pt, src_pt, out_pt;
 
   struct srcPt_t
   {
      double x[3];
      unsigned int index, proc;
   };
 
   struct outPt_t
   {
      double r[2], dist2;
      unsigned int index, code, el, proc;
      unsigned int mfem_el, loc_id, geom; // for triangles meshes.
   };
 
   {
      // Sync from device to host to send points to other ranks
      auto h_gsl_code = gsl_code.HostReadWrite(); // read now, write later
      auto h_pp = point_pos.HostRead();
 
      int index;
      struct srcPt_t *pt;
 
      array_init(struct srcPt_t, &hash_pt, points_cnt);
      pt = (struct srcPt_t *)hash_pt.ptr;
 
      Vector x(spacedim);
      for (index=0; index<points_cnt; ++index)
      {
         if (h_gsl_code[index] != CODE_INTERNAL)
         {
            for (int d=0; d<spacedim; ++d)
            {
               int idx = point_pos_ordering == 0 ?
                         index + d*points_cnt :
                         index*spacedim + d;
               x[d] = h_pp[idx];
            }
            const auto hi = hash_index_nd(DEV.gh_min, DEV.gh_fac, DEV.gh_nx, x);
            for (int d=0; d<spacedim; ++d)
            {
               pt->x[d] = x[d];
            }
            pt->index = index;
            pt->proc = hi % np;
            ++pt;
         }
      }
      hash_pt.n = pt - (struct srcPt_t *)hash_pt.ptr;
      sarray_transfer(struct srcPt_t, &hash_pt, proc, 1, DEV.cr);
   }
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
 
   /* look up points in hash cells, route to possible procs */
   {
      const unsigned int *const hash_offset = DEV.gh_offset.GetData();
      int count = 0;
      unsigned int *proc, *proc_p;
      const struct srcPt_t *p = (struct srcPt_t *)hash_pt.ptr,
                            *const pe = p + hash_pt.n;
      struct srcPt_t *q;
      Vector x(spacedim);
 
      for (; p!=pe; ++p)
      {
         for (int d = 0; d < spacedim; d++) { x[d] = p->x[d]; }
         const int hi = hash_index_nd(DEV.gh_min, DEV.gh_fac, DEV.gh_nx, x)/np;
         const int i = hash_offset[hi], ie = hash_offset[hi + 1];
         count += ie - i;
      }
 
      Array<unsigned int> proc_array(count);
      proc = proc_array.GetData();
      proc_p = proc;
      array_init(struct srcPt_t, &src_pt, count);
      q = (struct srcPt_t *)src_pt.ptr;
      p = (struct srcPt_t *)hash_pt.ptr;
      for (; p!=pe; ++p)
      {
         for (int d = 0; d < spacedim; d++) { x[d] = p->x[d]; }
         const int hi = hash_index_nd(DEV.gh_min, DEV.gh_fac, DEV.gh_nx, x)/np;
         int i        = hash_offset[hi];
         const int ie = hash_offset[hi + 1];
         for (; i!=ie; ++i)
         {
            const unsigned int pp = hash_offset[i];
            if (pp == p->proc)
            {
               continue;   /* don't send back to source proc */
            }
            *proc_p++ = pp;
            *q++      = *p;
         }
      }
 
      array_free(&hash_pt);
      src_pt.n = proc_p - proc;
 
      sarray_transfer_ext(struct srcPt_t, &src_pt, proc, sizeof(uint), DEV.cr);
   }
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
 
   /* look for other procs' points, send back */
   {
      int n = src_pt.n;
      const struct srcPt_t *spt;
      struct outPt_t *opt;
      array_init(struct outPt_t, &out_pt, n);
      out_pt.n = n;
      spt = (struct srcPt_t *)src_pt.ptr;
      opt = (struct outPt_t *)out_pt.ptr;
      for (; n; --n, ++spt, ++opt)
      {
         opt->index = spt->index;
         opt->proc  = spt->proc;
      }
      spt = (struct srcPt_t *)src_pt.ptr;
      opt = (struct outPt_t *)out_pt.ptr;
 
      n = out_pt.n;
      Vector gsl_ref_l(n*dim), gsl_dist_l(n);
      gsl_ref_l.UseDevice(use_dev);
      gsl_dist_l.UseDevice(use_dev);
 
      Vector point_pos_l(n*spacedim);
      point_pos_l.UseDevice(use_dev);
      auto pointl = point_pos_l.HostWrite();
 
      Array<unsigned int> gsl_code_l(n), gsl_elem_l(n);
 
      for (int point=0; point<n; ++point)
      {
         for (int d=0; d<spacedim; d++)
         {
            int idx = point_pos_ordering==0 ? point + d*n :
                      point*spacedim + d;
            pointl[idx] = spt[point].x[d];
         }
      }
 
      if (spacedim==2)
      {
         FindPointsEdgeLocal2(point_pos_l,
                              point_pos_ordering,
                              gsl_code_l,
                              gsl_elem_l,
                              gsl_ref_l,
                              gsl_dist_l,
                              n);
      }
      else
      {
         if (dim == 1)
         {
            FindPointsEdgeLocal3(point_pos_l,
                                 point_pos_ordering,
                                 gsl_code_l,
                                 gsl_elem_l,
                                 gsl_ref_l,
                                 gsl_dist_l,
                                 n);
 
         }
         else
         {
            FindPointsSurfLocal3(point_pos_l,
                                 point_pos_ordering,
                                 gsl_code_l,
                                 gsl_elem_l,
                                 gsl_ref_l,
                                 gsl_dist_l,
                                 n);
         }
      }
 
      gsl_ref_l.HostRead();
      gsl_dist_l.HostRead();
      gsl_code_l.HostRead();
      gsl_elem_l.HostRead();
      auto h_split_element_map = split_element_map.HostRead();
      auto h_split_element_index = split_element_index.HostRead();
      auto h_split_element_geom = split_element_geom.HostRead();
 
      // unpack arrays into opt
      for (int point=0; point<n; ++point)
      {
         opt[point].code = AsConst(gsl_code_l)[point];
         if (opt[point].code == CODE_NOT_FOUND)
         {
            continue;
         }
         opt[point].el   = AsConst(gsl_elem_l)[point];
         opt[point].dist2 = AsConst(gsl_dist_l)[point];
         for (int d = 0; d < dim; ++d)
         {
            opt[point].r[d] = AsConst(gsl_ref_l)[dim * point + d];
         }
         {
            const int loc_id = AsConst(gsl_elem_l)[point];
            opt[point].mfem_el = h_split_element_map[loc_id];
            opt[point].loc_id  = h_split_element_index[loc_id];
            opt[point].geom    = h_split_element_geom[loc_id];
         }
         // Note: we check if the point is on element border and mark it as
         // such. We do not mark points as CODE_INTERNAL because the found
         // solution could be interior to the element even when the point
         // is not on the edge/surface. This case is handled in the kernels.
         if (dim == 1)
         {
            real_t ipx = 0.5*opt[point].r[0]+0.5;
            if (ipx < rbtol || ipx > 1.0 - rbtol)
            {
               opt[point].code = CODE_BORDER;
            }
         }
         else if (dim == 2)
         {
            if (opt[point].geom == (int)Geometry::TRIANGLE)
            {
               // opt[point].el and r are for gslib.
               const int loc_id = opt[point].loc_id; // 0,1,2
               const double u = 0.5*opt[point].r[0] + 0.5;
               const double v = 0.5*opt[point].r[1] + 0.5;
               double tx, ty;
               MapSplitTriangleQuadToTriangle(loc_id, u, v, tx, ty);
               // Triangle edge or not.
               if (tx < rbtol || ty < rbtol || tx + ty > 1.0 - rbtol)
               {
                  opt[point].code = CODE_BORDER;
               }
            }
            else // quad
            {
               real_t ipx = 0.5*opt[point].r[0]+0.5;
               real_t ipy = 0.5*opt[point].r[1]+0.5;
               if (ipx < rbtol || ipx > 1.0 - rbtol ||
                   ipy < rbtol || ipy > 1.0 - rbtol)
               {
                  opt[point].code = CODE_BORDER;
               }
            }
         }
         if (opt[point].code == CODE_BORDER && opt[point].dist2 > bdr_tol)
         {
            opt[point].code = CODE_NOT_FOUND;
         }
      }
      array_free(&src_pt);
 
      /* group by code to eliminate unfound points */
      sarray_sort(struct outPt_t, opt, out_pt.n, code, 0, &DEV.cr->data);
 
      n = out_pt.n;
      while (n && opt[n-1].code == CODE_NOT_FOUND)
      {
         --n;
      }
      out_pt.n = n;
 
      sarray_transfer(struct outPt_t, &out_pt, proc, 1, DEV.cr);
   }
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
 
   auto h_gsl_code = gsl_code.HostReadWrite();
   auto h_gsl_dist = gsl_dist.HostReadWrite();
   auto h_gsl_ref = gsl_ref.HostReadWrite();
   auto h_gsl_mfem_ref = gsl_mfem_ref.HostReadWrite();
   auto h_gsl_elem = gsl_elem.HostReadWrite();
   auto h_gsl_mfem_elem = gsl_mfem_elem.HostReadWrite();
   auto h_gsl_proc = gsl_proc.HostReadWrite();
 
   // /* merge remote results with user data */
   {
      int n = out_pt.n;
      struct outPt_t *opt = (struct outPt_t *)out_pt.ptr;
      for (; n; --n, ++opt)
      {
         const int index = opt->index;
         if (h_gsl_code[index] == CODE_INTERNAL)
         {
            continue;
         }
         if ( h_gsl_code[index]==CODE_NOT_FOUND
              || opt->code==CODE_INTERNAL
              || opt->dist2<h_gsl_dist[index] )
         {
            for (int d=0; d<dim; ++d)
            {
               h_gsl_ref[dim*index + d] = opt->r[d];
            }
            h_gsl_dist[index]    = opt->dist2;
            h_gsl_proc[index]    = opt->proc;
            h_gsl_elem[index]    = opt->el;
            h_gsl_code[index]    = opt->code;
            h_gsl_mfem_elem[index] = opt->mfem_el;
 
            if (dim == 2 && opt->geom == (int)Geometry::TRIANGLE)
            {
               const int loc_id = opt->loc_id;
               const double u = 0.5*opt->r[0] + 0.5;
               const double v = 0.5*opt->r[1] + 0.5;
               MapSplitTriangleQuadToTriangle(loc_id, u, v,
                                              h_gsl_mfem_ref[dim*index + 0],
                                              h_gsl_mfem_ref[dim*index + 1]);
            }
            else
            {
               for (int d=0; d<dim; ++d)
               {
                  h_gsl_mfem_ref[dim*index + d] = 0.5*(opt->r[d] + 1.);
               }
            }
         }
      }
      array_free(&out_pt);
   }
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
}
 
void FindPointsGSLIB::InterpolateSurfBase(const Vector &field_in,
                                          Vector &field_out,
                                          const int nel,
                                          const int ncomp,
                                          const int dof1Dsol,
                                          const int field_out_ordering)
{
   bool use_dev = field_in.UseDevice();
   struct gslib::array src, outpt;
   int nlocal = 0;
   /* weed out unfound points, send out */
   Array<int> gsl_elem_temp;
   Vector gsl_ref_temp;
   Array<int> index_temp;
   {
      int index;
      const unsigned int *h_code = gsl_code.HostRead(),
                          *h_proc = gsl_proc.HostRead(),
                           *h_el   = gsl_elem.HostRead();
      const double *h_r = gsl_ref.HostRead();
 
      int numSend = 0;
 
      for (index=0; index<points_cnt; ++index)
      {
         numSend += (h_code[index] != CODE_NOT_FOUND &&
                     h_proc[index] != gsl_comm->id);
         nlocal += (h_code[index] != CODE_NOT_FOUND &&
                    h_proc[index] == gsl_comm->id);
      }
 
      gsl_elem_temp.SetSize(nlocal);
      gsl_elem_temp.HostWrite();
 
      gsl_ref_temp.SetSize(nlocal*dim);
      gsl_ref_temp.UseDevice(use_dev);
      gsl_ref_temp.HostWrite();
 
      index_temp.SetSize(nlocal);
 
      evalSrcPt_t *pt;
      array_init(evalSrcPt_t, &src, numSend);
      pt = (evalSrcPt_t *)src.ptr;
 
      int ctr = 0;
      for (index=0; index<points_cnt; ++index)
      {
         if (h_code[index]!=CODE_NOT_FOUND && h_proc[index]!=gsl_comm->id)
         {
            for (int d=0; d<dim; ++d)
            {
               pt->r[d] = h_r[index*dim + d];
            }
            pt->index = index;
            pt->proc = h_proc[index];
            pt->el = h_el[index];
            ++pt;
         }
         else if (h_code[index]!=CODE_NOT_FOUND && h_proc[index]==gsl_comm->id)
         {
            gsl_elem_temp[ctr] = h_el[index];
            for (int d=0; d<dim; ++d)
            {
               gsl_ref_temp(dim*ctr+d) = h_r[index*dim + d];
            }
            index_temp[ctr] = index;
            ctr++;
         }
      }
 
      src.n = pt - (evalSrcPt_t *)src.ptr;
      sarray_transfer(evalSrcPt_t, &src, proc, 1, cr);
   }
 
   //evaluate points that are already local
   {
      Vector interp_vals(nlocal*ncomp);
      interp_vals.UseDevice(use_dev);
 
      if (dim == 1)
      {
         InterpolateLocal1(field_in, gsl_elem_temp, gsl_ref_temp,
                           interp_vals, nlocal, ncomp, dof1Dsol);
      }
      else if (dim == 2)
      {
         InterpolateLocal2(field_in, gsl_elem_temp, gsl_ref_temp,
                           interp_vals, nlocal, ncomp, dof1Dsol);
 
      }
#ifdef MFEM_USE_MPI
      MPI_Barrier(gsl_comm->c);
#endif
 
      auto d_interp_vals = interp_vals.Read(use_dev);
      auto d_index_temp  = index_temp.Read(use_dev);
      auto d_field_out   = field_out.ReadWrite(use_dev); // no-op
      const int interp_offset = interp_vals.Size()/ncomp;
      InterpolateDeviceScatter(nlocal, d_index_temp, d_interp_vals,
                               d_field_out, interp_offset, ncomp, points_cnt,
                               field_out_ordering);
   }
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
 
   if (gsl_comm->np == 1)
   {
      array_free(&src);
      return;
   }
 
   // evaluate points locally
   {
      int n = src.n;
      const evalSrcPt_t *spt;
      evalOutPt_t *opt;
      spt = (evalSrcPt_t *)src.ptr;
 
      // Copy to host vector
      gsl_elem_temp.SetSize(n);
      gsl_elem_temp.HostWrite();
 
      gsl_ref_temp.SetSize(n*dim);
      gsl_ref_temp.HostWrite();
 
      spt = (evalSrcPt_t *)src.ptr;
      //   opt = (evalOutPt_t *)outpt.ptr;
      for (int i=0; i<n; i++, ++spt)
      {
         gsl_elem_temp[i] = spt->el;
         for (int d=0; d<dim; d++)
         {
            gsl_ref_temp(i*dim + d) = spt->r[d];
         }
      }
 
      Vector interp_vals(n*ncomp);
      interp_vals.UseDevice(use_dev);
      if (dim == 1)
      {
         InterpolateLocal1(field_in, gsl_elem_temp, gsl_ref_temp,
                           interp_vals, n, ncomp, dof1Dsol);
      }
      else if (dim == 2)
      {
         InterpolateLocal2(field_in, gsl_elem_temp, gsl_ref_temp,
                           interp_vals, n, ncomp, dof1Dsol);
      }
#ifdef MFEM_USE_MPI
      MPI_Barrier(gsl_comm->c);
#endif
      auto h_interp_vals = interp_vals.HostRead();
      auto h_field_out   = field_out.HostReadWrite();
 
      // Now the interpolated values need to be sent back component wise
      int Offset = interp_vals.Size()/ncomp;
      for (int i=0; i<ncomp; i++)
      {
         spt = (evalSrcPt_t *)src.ptr;
         array_init(evalOutPt_t, &outpt, n);
         outpt.n = n;
         opt = (evalOutPt_t *)outpt.ptr;
 
         for (int j=0; j<n; j++)
         {
            opt->index = spt->index;
            opt->proc = spt->proc;
            opt->out = h_interp_vals[j + Offset*i];
            spt++;
            opt++;
         }
 
         sarray_transfer(struct evalOutPt_t, &outpt, proc, 1, cr);
 
         opt = (evalOutPt_t *)outpt.ptr;
         for (size_t index = 0; index < outpt.n; index++)
         {
            int idx = field_out_ordering == Ordering::byNODES ?
                      opt->index + i*points_cnt :
                      opt->index*ncomp + i;
            h_field_out[idx] = opt->out;
            ++opt;
         }
         array_free(&outpt);
      }
      array_free(&src);
   }
   //finished evaluating points received from other processors.
}
 
#else
void FindPointsGSLIB::SetupDevice() {};
void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
                                         const int point_pos_ordering) {};
void FindPointsGSLIB::InterpolateOnDevice(const Vector &field_in_evec,
                                          Vector &field_out,
                                          const int nel, const int ncomp,
                                          const int dof1dsol,
                                          const int ordering) {};
 
void FindPointsGSLIB::FindPointsSurf(const Vector &point_pos,
                                     int point_pos_ordering)
{
   MFEM_ABORT("FindPointsGSLIB::FindPoints only supports surface meshes "
              "with GSLIB v1.0.9 or later.");
};
void FindPointsGSLIB::InterpolateSurfBase(const Vector &field_in_evec,
                                          Vector &field_out,
                                          const int nel, const int ncomp,
                                          const int dof1dsol,
                                          const int field_out_ordering) {};
#endif
 
void FindPointsGSLIB::FindPoints(Mesh &m, const Vector &point_pos,
                                 const int point_pos_ordering,
                                 const double bbox_rel_size_inc,
                                 const double newt_tol, const int npt_max)
{
   if (!setupflag || (mesh != &m))
   {
      Setup(m, bbox_rel_size_inc, newt_tol, npt_max);
   }
   FindPoints(point_pos, point_pos_ordering);
}
 
void FindPointsGSLIB::Interpolate(const Vector &point_pos,
                                  const GridFunction &field_in,
                                  Vector &field_out,
                                  const int point_pos_ordering)
{
   FindPoints(point_pos, point_pos_ordering);
   Interpolate(field_in, field_out);
}
 
void FindPointsGSLIB::Interpolate(const Vector &point_pos,
                                  const GridFunction &field_in,
                                  Vector &field_out,
                                  const int point_pos_ordering,
                                  const int field_out_ordering)
{
   FindPoints(point_pos, point_pos_ordering);
   Interpolate(field_in, field_out, field_out_ordering);
}
 
void FindPointsGSLIB::Interpolate(Mesh &m, const Vector &point_pos,
                                  const GridFunction &field_in,
                                  Vector &field_out,
                                  const int point_pos_ordering)
{
   FindPoints(m, point_pos, point_pos_ordering);
   Interpolate(field_in, field_out);
}
 
void FindPointsGSLIB::FreeData()
{
   FreeCrystal();
   if (!setupflag) { return; }
   if (dim == spacedim)
   {
#ifdef MFEM_USE_MPI
      if (!Mpi::IsFinalized())  // currently segfaults inside gslib otherwise
#endif
      {
         if (dim == 2)
         {
            findpts_free_2((gslib::findpts_data_2 *)this->fdataD);
         }
         else
         {
            findpts_free_3((gslib::findpts_data_3 *)this->fdataD);
         }
      }
   }
   gsl_code.DeleteAll();
   gsl_proc.DeleteAll();
   gsl_elem.DeleteAll();
   gsl_mesh.Destroy();
   gsl_ref.Destroy();
   gsl_dist.Destroy();
   for (int i = 0; i < 4; i++)
   {
      if (mesh_split[i]) { delete mesh_split[i]; mesh_split[i] = nullptr; }
      if (ir_split[i]) { delete ir_split[i]; ir_split[i] = nullptr; }
      if (ir_split_sol[i]) { delete ir_split_sol[i]; ir_split_sol[i] = nullptr; }
      if (fes_rst_map[i]) { delete fes_rst_map[i]; fes_rst_map[i] = nullptr; }
      if (gf_rst_map[i]) { delete gf_rst_map[i]; gf_rst_map[i] = nullptr; }
   }
   ir_split_sol_order = -1;
   if (fec_map_lin) { delete fec_map_lin; fec_map_lin = nullptr; }
   setupflag = false;
   DEV.setup_device = false;
   DEV.find_device  = false;
   points_cnt = -1;
}
 
void FindPointsGSLIB::SetupSplitMeshes()
{
   if (fec_map_lin == nullptr) { fec_map_lin = new H1_FECollection(1, dim); }
   if (mesh->Dimension() == 1)
   {
      mesh_split[0] = new Mesh(Mesh::MakeCartesian1D(1, Element::SEGMENT));
   }
   else if (mesh->Dimension() == 2)
   {
      int Nvert = 7;
      int NEsplit = 3;
      mesh_split[0] = new Mesh(2, Nvert, NEsplit, 0, 2);
 
      const double quad_v[7][2] =
      {
         {0, 0}, {0.5, 0}, {1, 0}, {0, 0.5},
         {1./3., 1./3.}, {0.5, 0.5}, {0, 1}
      };
      const int quad_e[3][4] =
      {
         {0, 1, 4, 3}, {1, 2, 5, 4}, {3, 4, 5, 6}
      };
 
      for (int j = 0; j < Nvert; j++)
      {
         mesh_split[0]->AddVertex(quad_v[j]);
      }
      for (int j = 0; j < NEsplit; j++)
      {
         int attribute = j + 1;
         mesh_split[0]->AddQuad(quad_e[j], attribute);
      }
      mesh_split[0]->FinalizeQuadMesh(1, 1, true);
 
      fes_rst_map[0] = new FiniteElementSpace(mesh_split[0], fec_map_lin, dim);
      gf_rst_map[0] = new GridFunction(fes_rst_map[0]);
      gf_rst_map[0]->UseDevice(false);
      const int npt = gf_rst_map[0]->Size()/dim;
      for (int k = 0; k < dim; k++)
      {
         for (int j = 0; j < npt; j++)
         {
            (*gf_rst_map[0])(j+k*npt) = quad_v[j][k];
         }
      }
 
      mesh_split[1] = new Mesh(Mesh::MakeCartesian2D(1, 1,
                                                     Element::QUADRILATERAL));
   }
   else if (dim == 3)
   {
      mesh_split[0] = new Mesh(Mesh::MakeCartesian3D(1, 1, 1,
                                                     Element::HEXAHEDRON));
      // Tetrahedron
      {
         int Nvert = 15;
         int NEsplit = 4;
         mesh_split[1] = new Mesh(3, Nvert, NEsplit, 0, 3);
 
         const double hex_v[15][3] =
         {
            {0, 0, 0.}, {1, 0., 0.}, {0., 1., 0.}, {0, 0., 1.},
            {0.5, 0., 0.}, {0.5, 0.5, 0.}, {0., 0.5, 0.},
            {0., 0., 0.5}, {0.5, 0., 0.5}, {0., 0.5, 0.5},
            {1./3., 0., 1./3.}, {1./3., 1./3., 1./3.}, {0, 1./3., 1./3.},
            {1./3., 1./3., 0}, {0.25, 0.25, 0.25}
         };
         const int hex_e[4][8] =
         {
            {7, 10, 4, 0, 12, 14, 13, 6},
            {10, 8, 1, 4, 14, 11, 5, 13},
            {14, 11, 5, 13, 12, 9, 2, 6},
            {7, 3, 8, 10, 12, 9, 11, 14}
         };
 
         for (int j = 0; j < Nvert; j++)
         {
            mesh_split[1]->AddVertex(hex_v[j]);
         }
         for (int j = 0; j < NEsplit; j++)
         {
            int attribute = j + 1;
            mesh_split[1]->AddHex(hex_e[j], attribute);
         }
         mesh_split[1]->FinalizeHexMesh(1, 1, true);
 
         fes_rst_map[1] = new FiniteElementSpace(mesh_split[1], fec_map_lin, dim);
         gf_rst_map[1] = new GridFunction(fes_rst_map[1]);
         gf_rst_map[1]->UseDevice(false);
         const int npt = gf_rst_map[1]->Size()/dim;
         for (int k = 0; k < dim; k++)
         {
            for (int j = 0; j < npt; j++)
            {
               (*gf_rst_map[1])(j+k*npt) = hex_v[j][k];
            }
         }
      }
      // Prism
      {
         int Nvert = 14;
         int NEsplit = 3;
         mesh_split[2] = new Mesh(3, Nvert, NEsplit, 0, 3);
 
         const double hex_v[14][3] =
         {
            {0, 0, 0}, {0.5, 0, 0}, {1, 0, 0}, {0, 0.5, 0},
            {1./3., 1./3., 0}, {0.5, 0.5, 0}, {0, 1, 0},
            {0, 0, 1}, {0.5, 0, 1}, {1, 0, 1}, {0, 0.5, 1},
            {1./3., 1./3., 1}, {0.5, 0.5, 1}, {0, 1, 1}
         };
         const int hex_e[3][8] =
         {
            {0, 1, 4, 3, 7, 8, 11, 10},
            {1, 2, 5, 4, 8, 9, 12, 11},
            {3, 4, 5, 6, 10, 11, 12, 13}
         };
 
         for (int j = 0; j < Nvert; j++)
         {
            mesh_split[2]->AddVertex(hex_v[j]);
         }
         for (int j = 0; j < NEsplit; j++)
         {
            int attribute = j + 1;
            mesh_split[2]->AddHex(hex_e[j], attribute);
         }
         mesh_split[2]->FinalizeHexMesh(1, 1, true);
 
         fes_rst_map[2] = new FiniteElementSpace(mesh_split[2], fec_map_lin, dim);
         gf_rst_map[2] = new GridFunction(fes_rst_map[2]);
         gf_rst_map[2]->UseDevice(false);
         const int npt = gf_rst_map[2]->Size()/dim;
         for (int k = 0; k < dim; k++)
         {
            for (int j = 0; j < npt; j++)
            {
               (*gf_rst_map[2])(j+k*npt) = hex_v[j][k];
            }
         }
      }
      // Pyramid
      {
         int Nvert = 23;
         int NEsplit = 8;
         mesh_split[3] = new Mesh(3, Nvert, NEsplit, 0, 3);
 
         const double hex_v[23][3] =
         {
            {0.0000, 0.0000, 0.0000}, {0.5000, 0.0000, 0.0000},
            {0.0000, 0.0000, 0.5000}, {0.3333, 0.0000, 0.3333},
            {0.0000, 0.5000, 0.0000}, {0.3333, 0.3333, 0.0000},
            {0.0000, 0.3333, 0.3333}, {0.2500, 0.2500, 0.2500},
            {1.0000, 0.0000, 0.0000}, {0.5000, 0.0000, 0.5000},
            {0.5000, 0.5000, 0.0000}, {0.3333, 0.3333, 0.3333},
            {0.0000, 1.0000, 0.0000}, {0.0000, 0.5000, 0.5000},
            {0.0000, 0.0000, 1.0000}, {1.0000, 0.5000, 0.0000},
            {0.6667, 0.3333, 0.3333}, {0.6667, 0.6667, 0.0000},
            {0.5000, 0.5000, 0.2500}, {1.0000, 1.0000, 0.0000},
            {0.5000, 0.5000, 0.5000}, {0.5000, 1.0000, 0.0000},
            {0.3333, 0.6667, 0.3333}
         };
         const int hex_e[8][8] =
         {
            {2, 3, 1, 0, 6, 7, 5, 4}, {3, 9, 8, 1, 7, 11, 10, 5},
            {7, 11, 10, 5, 6, 13, 12, 4}, {2, 14, 9, 3, 6, 13, 11, 7},
            {9, 16, 15, 8, 11, 18, 17, 10}, {16, 20, 19, 15, 18, 22, 21, 17},
            {18, 22, 21, 17, 11, 13, 12, 10}, {9, 14, 20, 16, 11, 13, 22, 18}
         };
 
         for (int j = 0; j < Nvert; j++)
         {
            mesh_split[3]->AddVertex(hex_v[j]);
         }
         for (int j = 0; j < NEsplit; j++)
         {
            int attribute = j + 1;
            mesh_split[3]->AddHex(hex_e[j], attribute);
         }
         mesh_split[3]->FinalizeHexMesh(1, 1, true);
 
         fes_rst_map[3] = new FiniteElementSpace(mesh_split[3], fec_map_lin, dim);
         gf_rst_map[3] = new GridFunction(fes_rst_map[3]);
         gf_rst_map[3]->UseDevice(false);
         const int npt = gf_rst_map[3]->Size()/dim;
         for (int k = 0; k < dim; k++)
         {
            for (int j = 0; j < npt; j++)
            {
               (*gf_rst_map[3])(j+k*npt) = hex_v[j][k];
            }
         }
      }
   }
}
 
void FindPointsGSLIB::SetupIntegrationRuleForSplitMesh(Mesh *meshin,
                                                       IntegrationRule *irule,
                                                       int order)
{
   H1_FECollection fec(order, dim);
   FiniteElementSpace nodal_fes(meshin, &fec, dim);
   GridFunction nodesplit(&nodal_fes);
   meshin->GetNodes(nodesplit);
   const int NEsplit = meshin->GetNE();
 
   const int dof_cnt = nodal_fes.GetTypicalFE()->GetDof(),
             pts_cnt = NEsplit * dof_cnt;
   Vector irlist(dim * pts_cnt);
 
   const TensorBasisElement *tbe =
      dynamic_cast<const TensorBasisElement *>(nodal_fes.GetTypicalFE());
   MFEM_VERIFY(tbe != nullptr, "TensorBasis FiniteElement expected.");
   const Array<int> &dof_map = tbe->GetDofMap();
 
   DenseMatrix pos(dof_cnt, dim);
   Vector posV(pos.Data(), dof_cnt * dim);
   Array<int> xdofs(dof_cnt * dim);
 
   // Create an IntegrationRule on the nodes of the reference submesh.
   MFEM_ASSERT(irule->GetNPoints() == pts_cnt, "IntegrationRule does not have"
               "the correct number of points.");
   int pt_id = 0;
   for (int i = 0; i < NEsplit; i++)
   {
      nodal_fes.GetElementVDofs(i, xdofs);
      nodesplit.GetSubVector(xdofs, posV);
      for (int j = 0; j < dof_cnt; j++)
      {
         for (int d = 0; d < dim; d++)
         {
            irlist(pts_cnt * d + pt_id) = pos(dof_map[j], d);
         }
         irule->IntPoint(pt_id).x = irlist(pt_id);
         if (dim >= 2)
         {
            irule->IntPoint(pt_id).y = irlist(pts_cnt + pt_id);
         }
         if (dim == 3)
         {
            irule->IntPoint(pt_id).z = irlist(2*pts_cnt + pt_id);
         }
         pt_id++;
      }
   }
}
 
void FindPointsGSLIB::SetupIntegrationRules(const int order,
                                            Array<IntegrationRule *> &ir_out)
{
   MFEM_VERIFY(mesh, "Setup FindPointsGSLIB with mesh first.");
   const int dof1D = order+1;
   dim = mesh->Dimension();
 
   if (dim == 1)
   {
      if (ir_out[0]) { delete ir_out[0]; ir_out[0] = nullptr; }
      ir_out[0] = new IntegrationRule(pow(dof1D, dim));
      SetupIntegrationRuleForSplitMesh(mesh_split[0], ir_out[0], order);
   }
   else if (dim == 2)
   {
      if (ir_out[0]) { delete ir_out[0]; ir_out[0] = nullptr; }
      ir_out[0] = new IntegrationRule(3*pow(dof1D, dim));
      SetupIntegrationRuleForSplitMesh(mesh_split[0], ir_out[0], order);
 
      if (ir_out[1]) { delete ir_out[1]; ir_out[1] = nullptr; }
      ir_out[1] = new IntegrationRule(pow(dof1D, dim));
      SetupIntegrationRuleForSplitMesh(mesh_split[1], ir_out[1], order);
   }
   else if (dim == 3)
   {
      if (ir_out[0]) { delete ir_out[0]; ir_out[0] = nullptr; }
      ir_out[0] = new IntegrationRule(pow(dof1D, dim));
      SetupIntegrationRuleForSplitMesh(mesh_split[0], ir_out[0], order);
 
      if (ir_out[1]) { delete ir_out[1]; ir_out[1] = nullptr; }
      ir_out[1] = new IntegrationRule(4*pow(dof1D, dim));
      SetupIntegrationRuleForSplitMesh(mesh_split[1], ir_out[1], order);
 
      if (ir_out[2]) { delete ir_out[2]; ir_out[2] = nullptr; }
      ir_out[2] = new IntegrationRule(3*pow(dof1D, dim));
      SetupIntegrationRuleForSplitMesh(mesh_split[2], ir_out[2], order);
 
      if (ir_out[3]) { delete ir_out[3]; ir_out[3] = nullptr; }
      ir_out[3] = new IntegrationRule(8*pow(dof1D, dim));
      SetupIntegrationRuleForSplitMesh(mesh_split[3], ir_out[3], order);
   }
}
 
void FindPointsGSLIB::SetupSplitMeshesAndIntegrationRules(const int order)
{
   MFEM_VERIFY(mesh, "Setup FindPointsGSLIB with mesh first.");
 
   SetupSplitMeshes();
   SetupIntegrationRules(order, ir_split);
 
   // Setup map for non tensor-product elements
   NE_split_total = 0;
   split_element_map.SetSize(0);
   split_element_index.SetSize(0);
   split_element_geom.SetSize(0);
   int NEsplit = 0;
   for (int e = 0; e < mesh->GetNE(); e++)
   {
      const Geometry::Type gt   = mesh->GetElement(e)->GetGeometryType();
      if (gt == Geometry::TRIANGLE || gt == Geometry::PRISM)
      {
         NEsplit = 3;
      }
      else if (gt == Geometry::TETRAHEDRON)
      {
         NEsplit = 4;
      }
      else if (gt == Geometry::PYRAMID)
      {
         NEsplit = 8;
      }
 
      else if (gt == Geometry::SEGMENT  ||
               gt == Geometry::SQUARE || gt == Geometry::CUBE)
      {
         NEsplit = 1;
      }
      else
      {
         MFEM_ABORT("Unsupported geometry type.");
      }
      NE_split_total += NEsplit;
      for (int i = 0; i < NEsplit; i++)
      {
         split_element_map.Append(e);
         split_element_index.Append(i);
         split_element_geom.Append((int)gt);
      }
   }
}
 
void FindPointsGSLIB::GetNodalValues(const GridFunction *gf_in,
                                     Vector &node_vals,
                                     const Array<IntegrationRule *> *ir_in,
                                     bool by_element) const
{
   const GridFunction *nodes     = gf_in;
   const FiniteElementSpace *fes = nodes->FESpace();
   const int NE                  = mesh->GetNE();
   const int vdim                = fes->GetVDim();
 
   const Array<IntegrationRule *> &ir = ir_in ? *ir_in : ir_split;
   IntegrationRule *ir_split_temp = nullptr;
 
   const int maxOrder = fes->GetMaxElementOrder();
   const int dof_1D =  maxOrder+1;
   const int pts_el = std::pow(dof_1D, dim);
   const int pts_cnt = NE_split_total * pts_el;
   node_vals.SetSize(vdim * pts_cnt);
 
   if (node_vals.UseDevice())
   {
      node_vals.HostWrite();
   }
 
   int gsl_mesh_pt_index = 0;
 
   for (int e = 0; e < NE; e++)
   {
      const FiniteElement *fe   = fes->GetFE(e);
      const Geometry::Type gt   = fe->GetGeomType();
      bool el_to_split = true;
      if (gt == Geometry::TRIANGLE)
      {
         ir_split_temp = ir[0];
      }
      else if (gt == Geometry::TETRAHEDRON)
      {
         ir_split_temp = ir[1];
      }
      else if (gt == Geometry::PRISM)
      {
         ir_split_temp = ir[2];
      }
      else if (gt == Geometry::PYRAMID)
      {
         ir_split_temp = ir[3];
      }
      else if (gt == Geometry::SQUARE)
      {
         ir_split_temp = ir[1];
         // "split" if input mesh is not a tensor basis or has mixed order
         el_to_split =
            gf_in->FESpace()->IsVariableOrder() ||
            dynamic_cast<const TensorBasisElement *>(fes->GetFE(e)) == nullptr;
      }
      else if (gt == Geometry::SEGMENT || gt == Geometry::CUBE)
      {
         ir_split_temp = ir[0];
         // "split" if input mesh is not a tensor basis or has mixed order
         el_to_split =
            gf_in->FESpace()->IsVariableOrder() ||
            dynamic_cast<const TensorBasisElement *>(fes->GetFE(e)) == nullptr;
      }
      else
      {
         MFEM_ABORT("Unsupported geometry type.");
      }
 
      if (el_to_split) // Triangle/Tet/Prism or Quads/Hex but variable order
      {
         // Fill node_vals with field values at split points.
         Vector locval(vdim);
         for (int i = 0; i < ir_split_temp->GetNPoints(); i++)
         {
            const IntegrationPoint &ip = ir_split_temp->IntPoint(i);
            nodes->GetVectorValue(e, ip, locval);
            if (by_element)
            {
               const int ei = gsl_mesh_pt_index / pts_el;
               const int pi  = gsl_mesh_pt_index % pts_el;
               for (int d = 0; d < vdim; d++)
               {
                  node_vals(ei * pts_el * vdim + d * pts_el + pi) = locval(d);
               }
            }
            else
            {
               for (int d = 0; d < vdim; d++)
               {
                  node_vals(pts_cnt*d + gsl_mesh_pt_index) = locval(d);
               }
            }
            gsl_mesh_pt_index++;
         }
      }
      else // Quad/Hex and constant polynomial order
      {
         const int dof_cnt_split = fe->GetDof();
 
         const TensorBasisElement *tbe =
            dynamic_cast<const TensorBasisElement *>(fes->GetFE(e));
         MFEM_VERIFY(tbe != nullptr, "TensorBasis FiniteElement expected.");
         Array<int> dof_map(dof_cnt_split);
         const Array<int> &dm = tbe->GetDofMap();
         if (dm.Size() > 0) { dof_map = dm; }
         else { for (int i = 0; i < dof_cnt_split; i++) { dof_map[i] = i; } }
 
         DenseMatrix pos(dof_cnt_split, vdim);
         Vector posV(pos.Data(), dof_cnt_split * vdim);
         Array<int> xdofs(dof_cnt_split * vdim);
 
         fes->GetElementVDofs(e, xdofs);
         nodes->GetSubVector(xdofs, posV);
         for (int j = 0; j < dof_cnt_split; j++)
         {
            if (by_element)
            {
               const int ei = gsl_mesh_pt_index / pts_el;
               const int pi  = gsl_mesh_pt_index % pts_el;
               for (int d = 0; d < vdim; d++)
               {
                  node_vals(ei * pts_el * vdim + d * pts_el + pi) =
                     pos(dof_map[j], d);
               }
            }
            else
            {
               for (int d = 0; d < vdim; d++)
               {
                  node_vals(pts_cnt * d + gsl_mesh_pt_index) = pos(dof_map[j], d);
               }
            }
            gsl_mesh_pt_index++;
         }
      }
   }
}
 
void FindPointsGSLIB::MapRefPosAndElemIndices()
{
   gsl_mfem_ref.SetSize(points_cnt*dim);
   gsl_mfem_elem.SetSize(points_cnt);
   gsl_mfem_ref = gsl_ref.HostRead();
   gsl_mfem_elem = gsl_elem;
 
   gsl_mfem_ref += 1.;  // map  [-1, 1] to [0, 2] to [0, 1]
   gsl_mfem_ref *= 0.5;
 
   int nptorig = points_cnt,
       npt = points_cnt;
 
   // Tolerance for point to be marked as on element edge/face based on the
   // obtained reference-space coordinates.
   double rbtol = 1e-12;
 
   GridFunction *gf_rst_map_temp = nullptr;
   int nptsend = 0;
 
   for (int index = 0; index < npt; index++)
   {
      if (gsl_code[index] != 2 && gsl_proc[index] != gsl_comm->id)
      {
         nptsend +=1;
      }
   }
 
   // Pack data to send via crystal router
   struct gslib::array *outpt = new gslib::array;
   struct out_pt { double r[3]; uint index, el, proc, code; };
   struct out_pt *pt;
   array_init(struct out_pt, outpt, nptsend);
   outpt->n=nptsend;
   pt = (struct out_pt *)outpt->ptr;
   for (int index = 0; index < npt; index++)
   {
      if (gsl_code[index] == 2 || gsl_proc[index] == gsl_comm->id)
      {
         continue;
      }
      for (int d = 0; d < dim; ++d)
      {
         pt->r[d]= gsl_mfem_ref(index*dim + d);
      }
      pt->index = index;
      pt->proc  = gsl_proc[index];
      pt->el    = gsl_elem[index];
      pt->code  = gsl_code[index];
      ++pt;
   }
 
   // Transfer data to target MPI ranks
   sarray_transfer(struct out_pt, outpt, proc, 1, cr);
 
   // Map received points
   npt = outpt->n;
   pt = (struct out_pt *)outpt->ptr;
   for (int index = 0; index < npt; index++)
   {
      IntegrationPoint ip;
      ip.Set3(&pt->r[0]);
      const int elem = pt->el;
      const int mesh_elem = split_element_map[elem];
      const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(mesh_elem);
 
      const Geometry::Type gt = fe->GetGeomType();
      pt->el = mesh_elem;
 
      if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
      {
         // check if it is on element boundary
         pt->code = Geometry::CheckPoint(gt, ip, -rbtol) ? 0 : 1;
         ++pt;
         continue;
      }
      else if (gt == Geometry::TRIANGLE)
      {
         gf_rst_map_temp = gf_rst_map[0];
      }
      else if (gt == Geometry::TETRAHEDRON)
      {
         gf_rst_map_temp = gf_rst_map[1];
      }
      else if (gt == Geometry::PRISM)
      {
         gf_rst_map_temp = gf_rst_map[2];
      }
      else if (gt == Geometry::PYRAMID)
      {
         gf_rst_map_temp = gf_rst_map[3];
      }
 
      int local_elem = split_element_index[elem];
      Vector mfem_ref(dim);
      // map to rst of macro element
      gf_rst_map_temp->GetVectorValue(local_elem, ip, mfem_ref);
 
      for (int d = 0; d < dim; d++)
      {
         pt->r[d] = mfem_ref(d);
      }
 
      // check if point is on element boundary
      ip.Set3(&pt->r[0]);
      pt->code = Geometry::CheckPoint(gt, ip, -rbtol) ? 0 : 1;
      ++pt;
   }
 
   // Transfer data back to source MPI rank
   sarray_transfer(struct out_pt, outpt, proc, 1, cr);
   npt = outpt->n;
 
   // First copy mapped information for points on other procs
   pt = (struct out_pt *)outpt->ptr;
   for (int index = 0; index < npt; index++)
   {
      gsl_mfem_elem[pt->index] = pt->el;
      for (int d = 0; d < dim; d++)
      {
         gsl_mfem_ref(d + pt->index*dim) = pt->r[d];
      }
      gsl_code[pt->index] = pt->code;
      ++pt;
   }
   array_free(outpt);
   delete outpt;
 
   // Now map information for points on the same proc
   for (int index = 0; index < nptorig; index++)
   {
      if (gsl_code[index] != 2 && gsl_proc[index] == gsl_comm->id)
      {
 
         IntegrationPoint ip;
         Vector mfem_ref(gsl_mfem_ref.GetData()+index*dim, dim);
         ip.Set2(mfem_ref.GetData());
         if (dim == 3) { ip.z = mfem_ref(2); }
 
         const int elem = gsl_elem[index];
         const int mesh_elem = split_element_map[elem];
         const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(mesh_elem);
         const Geometry::Type gt = fe->GetGeomType();
         gsl_mfem_elem[index] = mesh_elem;
         if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
         {
            gsl_code[index] = Geometry::CheckPoint(gt, ip, -rbtol) ? 0 : 1;
            continue;
         }
         else if (gt == Geometry::TRIANGLE)
         {
            gf_rst_map_temp = gf_rst_map[0];
         }
         else if (gt == Geometry::TETRAHEDRON)
         {
            gf_rst_map_temp = gf_rst_map[1];
         }
         else if (gt == Geometry::PRISM)
         {
            gf_rst_map_temp = gf_rst_map[2];
         }
         else if (gt == Geometry::PYRAMID)
         {
            gf_rst_map_temp = gf_rst_map[3];
         }
 
         int local_elem = split_element_index[elem];
         gf_rst_map_temp->GetVectorValue(local_elem, ip, mfem_ref);
 
         // Check if the point is on element boundary
         ip.Set2(mfem_ref.GetData());
         if (dim == 3) { ip.z = mfem_ref(2); }
         gsl_code[index]  = Geometry::CheckPoint(gt, ip, -rbtol) ? 0 : 1;
      }
   }
}
 
void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
                                  Vector &field_out)
{
   Interpolate(field_in, field_out, field_in.FESpace()->GetOrdering());
}
 
void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
                                  Vector &field_out,
                                  const int field_out_ordering)
{
   MFEM_VERIFY(setupflag, "FindPointsGSLIB::Setup must be called first.");
   if (dim != spacedim)
   {
      InterpolateSurf(field_in, field_out, field_out_ordering);
      return;
   }
   const int  gf_order   = field_in.FESpace()->GetMaxElementOrder(),
              mesh_order = mesh->GetNodalFESpace()->GetMaxElementOrder();
 
   const FiniteElementCollection *fec_in =  field_in.FESpace()->FEColl();
   const H1_FECollection *fec_h1 = dynamic_cast<const H1_FECollection *>(fec_in);
   const L2_FECollection *fec_l2 = dynamic_cast<const L2_FECollection *>(fec_in);
 
   bool tensor_product_only = mesh->GetNE() == 0 ||
                              (mesh->GetNumGeometries(dim) == 1 &&
                               (mesh->GetElementType(0)==Element::QUADRILATERAL ||
                                mesh->GetElementType(0) == Element::HEXAHEDRON));
#ifdef MFEM_USE_MPI
   MPI_Allreduce(MPI_IN_PLACE, &tensor_product_only, 1, MFEM_MPI_CXX_BOOL,
                 MPI_LAND, gsl_comm->c);
#endif
 
   bool field_in_on_dev = Device::IsEnabled() && field_in.UseDevice();
   bool field_out_on_dev = Device::IsEnabled() && field_out.UseDevice();
 
   if (field_in_on_dev && fec_h1 &&
       !field_in.FESpace()->IsVariableOrder() && tensor_product_only)
   {
#if GSLIB_RELEASE_VERSION == 10007
      if (!gpu_to_cpu_fallback)
      {
         MFEM_ABORT("Either update to gslib v1.0.9 for GPU support "
                    "or use SetGPUtoCPUFallback to use host-functions. See "
                    "INSTALL for instructions to update GSLIB");
      }
#else
      MFEM_VERIFY(fec_h1->GetBasisType() == BasisType::GaussLobatto,
                  "basis not supported");
      Vector node_vals;
      const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
      const Operator *R = field_in.FESpace()->GetElementRestriction(ordering);
      node_vals.UseDevice(true);
      node_vals.SetSize(R->Height(), Device::GetMemoryType());
      R->Mult(field_in, node_vals);
      // GetNodalValues(&field_in, node_vals);
 
      const int ncomp  = field_in.FESpace()->GetVDim();
      const int maxOrder = field_in.FESpace()->GetMaxElementOrder();
 
      InterpolateOnDevice(node_vals, field_out, NE_split_total, ncomp,
                          maxOrder+1, field_out_ordering);
      return;
#endif
   }
   // The remaining interpolation paths are host-only. They access field_out
   // through explicit host pointers and address device validity before
   // returning (if field_out or field_in is on device).
   field_in.HostRead();
 
   // Read GSLIB data to host in case it was done on device.
   auto h_gsl_code = gsl_code.HostRead();
   gsl_elem.HostRead();
   gsl_mfem_elem.HostRead();
   gsl_ref.HostRead();
   gsl_mfem_ref.HostRead();
   gsl_proc.HostRead();
 
   if (fec_h1 && gf_order == mesh_order &&
       fec_h1->GetBasisType() == BasisType::GaussLobatto &&
       field_in.FESpace()->IsVariableOrder() ==
       mesh->GetNodalFESpace()->IsVariableOrder())
   {
      InterpolateH1(field_in, field_out, field_out_ordering);
      if (field_in_on_dev || field_out_on_dev)
      {
         field_out.ReadWrite(); // also internally calls mem.UseDevice(true);
      }
      return;
   }
   else
   {
      InterpolateGeneral(field_in, field_out, field_out_ordering);
      if (!fec_l2 || avgtype == AvgType::NONE)
      {
         if (field_in_on_dev || field_out_on_dev)
         {
            field_out.ReadWrite(); // also internally calls mem.UseDevice(true);
         }
         return;
      }
   }
 
   // For points on element borders, project the L2 GridFunction to H1 and
   // re-interpolate.
   if (fec_l2)
   {
      Array<int> indl2;
      for (int i = 0; i < points_cnt; i++)
      {
         if (h_gsl_code[i] == 1) { indl2.Append(i); }
      }
      int borderPts = indl2.Size();
#ifdef MFEM_USE_MPI
      MPI_Allreduce(MPI_IN_PLACE, &borderPts, 1, MPI_INT, MPI_SUM, gsl_comm->c);
#endif
      if (borderPts == 0)
      {
         if (field_in_on_dev || field_out_on_dev)
         {
            field_out.ReadWrite(); // also internally calls mem.UseDevice(true);
         }
         return; // no points on element borders
      }
 
      Vector field_out_l2(field_out.Size());
      VectorGridFunctionCoefficient field_in_dg(&field_in);
      int gf_order_h1 = std::max(gf_order, 1); // H1 should be at least order 1
      H1_FECollection fec(gf_order_h1, dim);
      const int ncomp = field_in.FESpace()->GetVDim();
 
      std::unique_ptr<FiniteElementSpace> fes;
      std::unique_ptr<GridFunction> field_in_h1;
#ifdef MFEM_USE_MPI
      if (auto *pmesh = dynamic_cast<ParMesh*>(mesh))
      {
         fes = std::make_unique<ParFiniteElementSpace>(
                  pmesh, &fec, ncomp, field_in.FESpace()->GetOrdering());
 
         field_in_h1 = std::make_unique<ParGridFunction>(
                          static_cast<ParFiniteElementSpace*>(fes.get()));
      }
      else
#endif
      {
         fes = std::make_unique<FiniteElementSpace>(
                  mesh, &fec, ncomp, field_in.FESpace()->GetOrdering());
 
         field_in_h1 = std::make_unique<GridFunction>(fes.get());
      }
 
      field_in_h1->UseDevice(false);
 
      if (avgtype == AvgType::ARITHMETIC)
      {
         field_in_h1->ProjectDiscCoefficient(field_in_dg, GridFunction::ARITHMETIC);
      }
      else if (avgtype == AvgType::HARMONIC)
      {
         field_in_h1->ProjectDiscCoefficient(field_in_dg, GridFunction::HARMONIC);
      }
      else
      {
         MFEM_ABORT("Invalid averaging type.");
      }
 
      if (gf_order_h1 == mesh_order) // basis is GaussLobatto by default
      {
         InterpolateH1(*field_in_h1, field_out_l2, field_out_ordering);
      }
      else
      {
         InterpolateGeneral(*field_in_h1, field_out_l2, field_out_ordering);
      }
 
      // Copy interpolated values for the points on element border
      real_t *h_field_out = field_out.HostReadWrite();
      const real_t *h_field_out_l2 = field_out_l2.HostRead();
      for (int j = 0; j < ncomp; j++)
      {
         for (int i = 0; i < indl2.Size(); i++)
         {
            int idx = field_out_ordering == Ordering::byNODES?
                      indl2[i] + j*points_cnt:
                      indl2[i]*ncomp + j;
            h_field_out[idx] = h_field_out_l2[idx];
         }
      }
      if (field_in_on_dev || field_out_on_dev)
      {
         field_out.ReadWrite(); // also internally calls mem.UseDevice(true);
      }
   }
}
 
void FindPointsGSLIB::InterpolateSurf(const GridFunction &field_in,
                                      Vector &field_out)
{
   InterpolateSurf(field_in, field_out, field_in.FESpace()->GetOrdering());
}
 
void FindPointsGSLIB::InterpolateSurf(const GridFunction &field_in,
                                      Vector &field_out,
                                      const int field_out_ordering)
{
   MFEM_VERIFY(setupflag, "FindPointsGSLIB::Setup must be called first.");
#if GSLIB_RELEASE_VERSION < 10009
   MFEM_ABORT("Update to gslib v1.0.9 for surface mesh support.");
#endif
   const FiniteElementCollection *fec_in =  field_in.FESpace()->FEColl();
   const H1_FECollection *fec_h1 = dynamic_cast<const H1_FECollection *>
                                   (fec_in);
   MFEM_VERIFY(fec_h1,"Only h1 functions supported for surface meshes.");
   MFEM_VERIFY(fec_h1->GetBasisType() == BasisType::GaussLobatto,
               "basis not supported");
   MFEM_VERIFY(dim < spacedim, "InterpolateSurf is only for surface meshes.");
 
   bool has_split_elems = false;
   bool supported_surf_elem = true;
   for (int e = 0; e < mesh->GetNE(); e++)
   {
      const Element::Type t = mesh->GetElementType(e);
      if (t == Element::TRIANGLE) { has_split_elems = true; }
      else if (t != Element::SEGMENT && t != Element::QUADRILATERAL)
      {
         supported_surf_elem = false;
      }
   }
#ifdef MFEM_USE_MPI
   MPI_Allreduce(MPI_IN_PLACE, &supported_surf_elem, 1, MFEM_MPI_CXX_BOOL,
                 MPI_LAND, gsl_comm->c);
#endif
   MFEM_VERIFY(supported_surf_elem,
               "FindPointsGSLIB surface mesh support: only SEGMENT, "
               "QUADRILATERAL, and TRIANGLE elements are supported.");
   MFEM_VERIFY(dim < spacedim, "InterpolateSurf only supports surface meshes.");
 
   bool use_dev = field_in.UseDevice();
 
   if (!field_in.FESpace()->IsVariableOrder())
   {
      Vector node_vals;
 
      if (has_split_elems)
      {
         const int fieldOrder = field_in.FESpace()->GetMaxElementOrder();
         const int meshOrder  = DEV.dof1d - 1;
 
         const Array<IntegrationRule *> *ir_to_use;
         if (fieldOrder == meshOrder)
         {
            ir_to_use = &ir_split;
         }
         else
         {
            if (fieldOrder != ir_split_sol_order)
            {
               SetupIntegrationRules(fieldOrder, ir_split_sol);
               ir_split_sol_order = fieldOrder;
            }
            ir_to_use = &ir_split_sol;
         }
         node_vals.UseDevice(use_dev);
         GetNodalValues(&field_in, node_vals, ir_to_use, true);
      }
      else
      {
         node_vals.UseDevice(use_dev);
         const ElementDofOrdering ord = ElementDofOrdering::LEXICOGRAPHIC;
         const Operator *R            = field_in.FESpace()->GetElementRestriction(ord);
         node_vals.SetSize(R->Height());
         R->Mult(field_in, node_vals); // layout: [NEL][VDIM][N^D]
      }
 
      const int ncomp    = field_in.FESpace()->GetVDim();
      const int maxOrder = field_in.FESpace()->GetMaxElementOrder();
      DEV.dof1d_sol      =  maxOrder+1;
      DEV.gll1d_sol.UseDevice(use_dev);
      DEV.gll1d_sol.SetSize(DEV.dof1d_sol);
      DEV.lagcoeff_sol.UseDevice(use_dev);
      DEV.lagcoeff_sol.SetSize(DEV.dof1d_sol);
      if (DEV.dof1d_sol != DEV.dof1d)
      {
         auto h_gll1d_sol = DEV.gll1d_sol.HostWrite();
         auto h_lagcoeff_sol = DEV.lagcoeff_sol.HostWrite();
         gslib::lobatto_nodes(h_gll1d_sol, DEV.dof1d_sol);
         gslib::gll_lag_setup(h_lagcoeff_sol, DEV.dof1d_sol);
      }
      else
      {
         auto h_gll1d_sol = DEV.gll1d_sol.HostWrite();
         auto h_lagcoeff_sol = DEV.lagcoeff_sol.HostWrite();
         auto h_gll1d = DEV.gll1d.HostRead();
         auto h_lagcoeff = DEV.lagcoeff.HostRead();
         for (int i = 0; i < DEV.dof1d_sol; i++)
         {
            h_gll1d_sol[i] = h_gll1d[i];
            h_lagcoeff_sol[i] = h_lagcoeff[i];
         }
      }
 
      field_out.SetSize(points_cnt*ncomp);
      field_out.UseDevice(use_dev);
      field_out = default_interp_value;
 
      InterpolateSurfBase(node_vals, field_out, NE_split_total, ncomp,
                          DEV.dof1d_sol, field_out_ordering);
      return;
   }
   else
   {
      MFEM_ABORT("Variable order functions not supported for surface meshes!!");
   }
}
 
void FindPointsGSLIB::InterpolateH1(const GridFunction &field_in,
                                    Vector &field_out,
                                    const int field_out_ordering)
{
   FiniteElementSpace ind_fes(mesh, field_in.FESpace()->FEColl());
   if (field_in.FESpace()->IsVariableOrder())
   {
      for (int e = 0; e < ind_fes.GetMesh()->GetNE(); e++)
      {
         ind_fes.SetElementOrder(e, field_in.FESpace()->GetElementOrder(e));
      }
      ind_fes.Update(false);
   }
   GridFunction field_in_scalar(&ind_fes);
   field_in_scalar.UseDevice(false);
   Vector node_vals;
 
   const int ncomp      = field_in.FESpace()->GetVDim(),
             points_fld = field_in.Size() / ncomp;
   MFEM_VERIFY(points_cnt == gsl_code.Size(),
               "FindPointsGSLIB::InterpolateH1: Inconsistent size of gsl_code");
 
   field_out.SetSize(points_cnt*ncomp);
   real_t *h_field_out = field_out.HostWrite();
   std::fill(h_field_out, h_field_out + field_out.Size(),
             default_interp_value);
 
   for (int i = 0; i < ncomp; i++)
   {
      const int dataptrin  = i*points_fld,
                dataptrout = i*points_cnt;
      if (field_in.FESpace()->GetOrdering() == Ordering::byNODES)
      {
         field_in_scalar.NewDataAndSize(field_in.GetData()+dataptrin,
                                        points_fld);
      }
      else
      {
         for (int j = 0; j < points_fld; j++)
         {
            field_in_scalar(j) = field_in(i + j*ncomp);
         }
      }
      GetNodalValues(&field_in_scalar, node_vals);
 
      if (dim==2)
      {
         findpts_eval_2(h_field_out+dataptrout, sizeof(double),
                        gsl_code.GetData(),       sizeof(unsigned int),
                        gsl_proc.GetData(),    sizeof(unsigned int),
                        gsl_elem.GetData(),    sizeof(unsigned int),
                        gsl_ref.GetData(),     sizeof(double) * dim,
                        points_cnt, node_vals.GetData(),
                        (gslib::findpts_data_2 *)this->fdataD);
      }
      else
      {
         findpts_eval_3(h_field_out+dataptrout, sizeof(double),
                        gsl_code.GetData(),       sizeof(unsigned int),
                        gsl_proc.GetData(),    sizeof(unsigned int),
                        gsl_elem.GetData(),    sizeof(unsigned int),
                        gsl_ref.GetData(),     sizeof(double) * dim,
                        points_cnt, node_vals.GetData(),
                        (gslib::findpts_data_3 *)this->fdataD);
      }
   }
   if (field_out_ordering == Ordering::byVDIM)
   {
      Vector field_out_temp(field_out.Size()); // host-only
      field_out_temp = h_field_out;
      for (int i = 0; i < ncomp; i++)
      {
         for (int j = 0; j < points_cnt; j++)
         {
            h_field_out[i + j*ncomp] = field_out_temp(j + i*points_cnt);
         }
      }
   }
}
 
void FindPointsGSLIB::InterpolateGeneral(const GridFunction &field_in,
                                         Vector &field_out,
                                         const int field_out_ordering)
{
   int ncomp   = field_in.VectorDim(),
       nptorig = points_cnt,
       npt     = points_cnt;
 
   field_out.SetSize(points_cnt*ncomp);
   real_t *h_field_out = field_out.HostWrite();
   std::fill(h_field_out, h_field_out + field_out.Size(),
             default_interp_value);
 
   // Get host read pointers in case FindPoints was done on device.
   auto h_gsl_code = gsl_code.HostRead();
   auto h_gsl_mfem_ref = gsl_mfem_ref.HostRead();
   auto h_gsl_mfem_elem = gsl_mfem_elem.HostRead();
   auto h_gsl_proc = gsl_proc.HostRead();
 
   if (gsl_comm->np == 1) // serial
   {
      for (int index = 0; index < npt; index++)
      {
         if (h_gsl_code[index] == 2) { continue; }
         IntegrationPoint ip;
         ip.Set(h_gsl_mfem_ref + index*dim, dim);
         Vector localval(ncomp);
         field_in.GetVectorValue(h_gsl_mfem_elem[index], ip, localval);
         if (field_out_ordering == Ordering::byNODES)
         {
            for (int i = 0; i < ncomp; i++)
            {
               h_field_out[index + i*npt] = localval(i);
            }
         }
         else //byVDIM
         {
            for (int i = 0; i < ncomp; i++)
            {
               h_field_out[index*ncomp + i] = localval(i);
            }
         }
      }
   }
   else // parallel
   {
      // Determine number of points to be sent
      int nptsend = 0;
      for (int index = 0; index < npt; index++)
      {
         if (h_gsl_code[index] != 2) { nptsend +=1; }
      }
 
      // Pack data to send via crystal router
      struct gslib::array *outpt = new gslib::array;
      struct out_pt { double r[3], ival; uint index, el, proc; };
      struct out_pt *pt;
      array_init(struct out_pt, outpt, nptsend);
      outpt->n=nptsend;
      pt = (struct out_pt *)outpt->ptr;
      for (int index = 0; index < npt; index++)
      {
         if (h_gsl_code[index] == 2) { continue; }
         for (int d = 0; d < dim; ++d)
         {
            pt->r[d]= h_gsl_mfem_ref[index*dim + d];
         }
         pt->index = index;
         pt->proc  = h_gsl_proc[index];
         pt->el    = h_gsl_mfem_elem[index];
         ++pt;
      }
 
      // Transfer data to target MPI ranks
      sarray_transfer(struct out_pt, outpt, proc, 1, cr);
 
      if (ncomp == 1)
      {
         // Interpolate the grid function
         npt = outpt->n;
         pt = (struct out_pt *)outpt->ptr;
         for (int index = 0; index < npt; index++)
         {
            IntegrationPoint ip;
            ip.Set3(&pt->r[0]);
            pt->ival = field_in.GetValue(pt->el, ip, 1);
            ++pt;
         }
 
         // Transfer data back to source MPI rank
         sarray_transfer(struct out_pt, outpt, proc, 1, cr);
         npt = outpt->n;
         pt = (struct out_pt *)outpt->ptr;
         for (int index = 0; index < npt; index++)
         {
            h_field_out[pt->index] = pt->ival;
            ++pt;
         }
         array_free(outpt);
         delete outpt;
      }
      else // ncomp > 1
      {
         // Interpolate data and store in a Vector
         npt = outpt->n;
         pt = (struct out_pt *)outpt->ptr;
         Vector vec_int_vals(npt*ncomp);
         for (int index = 0; index < npt; index++)
         {
            IntegrationPoint ip;
            ip.Set3(&pt->r[0]);
            Vector localval(vec_int_vals.GetData()+index*ncomp, ncomp);
            field_in.GetVectorValue(pt->el, ip, localval);
            ++pt;
         }
 
         // Save index and proc data in a struct
         struct gslib::array *savpt = new gslib::array;
         struct sav_pt { uint index, proc; };
         struct sav_pt *spt;
         array_init(struct sav_pt, savpt, npt);
         savpt->n=npt;
         spt = (struct sav_pt *)savpt->ptr;
         pt  = (struct out_pt *)outpt->ptr;
         for (int index = 0; index < npt; index++)
         {
            spt->index = pt->index;
            spt->proc  = pt->proc;
            ++pt; ++spt;
         }
 
         array_free(outpt);
         delete outpt;
 
         // Copy data from save struct to send struct and send component wise
         struct gslib::array *sendpt = new gslib::array;
         struct send_pt { double ival; uint index, proc; };
         struct send_pt *sdpt;
         for (int j = 0; j < ncomp; j++)
         {
            array_init(struct send_pt, sendpt, npt);
            sendpt->n=npt;
            spt  = (struct sav_pt *)savpt->ptr;
            sdpt = (struct send_pt *)sendpt->ptr;
            for (int index = 0; index < npt; index++)
            {
               sdpt->index = spt->index;
               sdpt->proc  = spt->proc;
               sdpt->ival  = vec_int_vals(j + index*ncomp);
               ++sdpt; ++spt;
            }
 
            sarray_transfer(struct send_pt, sendpt, proc, 1, cr);
            sdpt = (struct send_pt *)sendpt->ptr;
            for (int index = 0; index < static_cast<int>(sendpt->n); index++)
            {
               int idx = field_out_ordering == Ordering::byNODES ?
                         sdpt->index + j*nptorig :
                         sdpt->index*ncomp + j;
               h_field_out[idx] = sdpt->ival;
               ++sdpt;
            }
            array_free(sendpt);
         }
         array_free(savpt);
         delete sendpt;
         delete savpt;
      } // ncomp > 1
   } // parallel
}
 
Array<unsigned int> FindPointsGSLIB::GetPointsNotFoundIndices() const
{
   Array<unsigned int> nf_idxs;
   auto h_gsl_code = gsl_code.HostRead();
   for (int i = 0; i < gsl_code.Size(); i++)
   {
      if (h_gsl_code[i] == 2)
      {
         nf_idxs.Append(i);
      }
   }
   return nf_idxs;
}
 
void FindPointsGSLIB::DistributePointInfoToOwningMPIRanks(
   Array<unsigned int> &recv_elem, Vector &recv_ref,
   Array<unsigned int> &recv_code)
{
   MFEM_VERIFY(points_cnt >= 0,
               "Invalid size. Please make sure to call FindPoints method "
               "before calling this function.");
 
   // Pack data to send via crystal router
   struct gslib::array *outpt = new gslib::array;
 
   struct out_pt { double rst[3]; uint index, elem, proc, code; };
   struct out_pt *pt;
   array_init(struct out_pt, outpt, points_cnt);
   outpt->n=points_cnt;
   pt = (struct out_pt *)outpt->ptr;
 
   for (int index = 0; index < points_cnt; index++)
   {
      pt->index = index;
      pt->elem = gsl_mfem_elem[index];
      pt->proc  = gsl_proc[index];
      pt->code = gsl_code[index];
      for (int d = 0; d < dim; ++d)
      {
         pt->rst[d]= gsl_mfem_ref(index*dim + d);
      }
      ++pt;
   }
 
   // Transfer data to target MPI ranks
   sarray_transfer(struct out_pt, outpt, proc, 1, cr);
 
   // Store received data
   const int points_recv = outpt->n;
   recv_proc.SetSize(points_recv);
   recv_elem.SetSize(points_recv);
   recv_index.SetSize(points_recv);
   recv_code.SetSize(points_recv);
   recv_ref.SetSize(points_recv*dim);
 
   pt = (struct out_pt *)outpt->ptr;
   for (int index = 0; index < points_recv; index++)
   {
      recv_index[index] = pt->index;
      recv_elem[index] = pt->elem;
      recv_proc[index] = pt->proc;
      recv_code[index] = pt->code;
      for (int d = 0; d < dim; ++d)
      {
         recv_ref(index*dim + d)= pt->rst[d];
      }
      ++pt;
   }
 
   array_free(outpt);
   delete outpt;
}
 
void FindPointsGSLIB::DistributeInterpolatedValues(const Vector &int_vals,
                                                   const int vdim,
                                                   const int ordering,
                                                   Vector &field_out) const
{
   const int points_recv = recv_index.Size();;
   MFEM_VERIFY(points_recv == 0 ||
               int_vals.Size() % points_recv == 0,
               "Incompatible size. Please return interpolated values"
               "corresponding to points received using"
               "SendCoordinatesToOwningProcessors.");
   field_out.SetSize(points_cnt*vdim);
 
   for (int v = 0; v < vdim; v++)
   {
      // Pack data to send via crystal router
      struct gslib::array *outpt = new gslib::array;
      struct out_pt { double val; uint index, proc; };
      struct out_pt *pt;
      array_init(struct out_pt, outpt, points_recv);
      outpt->n=points_recv;
      pt = (struct out_pt *)outpt->ptr;
      for (int index = 0; index < points_recv; index++)
      {
         pt->index = recv_index[index];
         pt->proc  = recv_proc[index];
         pt->val = ordering == Ordering::byNODES ?
                   int_vals(index + v*points_recv) :
                   int_vals(index*vdim + v);
         ++pt;
      }
 
      // Transfer data to target MPI ranks
      sarray_transfer(struct out_pt, outpt, proc, 1, cr);
 
      // Store received data
      MFEM_VERIFY(outpt->n == static_cast<size_t>(points_cnt),
                  "Incompatible size. Number of points "
                  "received does not match the number of points originally "
                  "found using FindPoints.");
 
      pt = (struct out_pt *)outpt->ptr;
      for (int index = 0; index < points_cnt; index++)
      {
         int idx = ordering == Ordering::byNODES ?
                   pt->index + v*points_cnt :
                   pt->index*vdim + v;
         field_out(idx) = pt->val;
         ++pt;
      }
 
      array_free(outpt);
      delete outpt;
   }
}
 
void FindPointsGSLIB::GetAxisAlignedBoundingBoxes(Vector &aabb) const
{
   MFEM_VERIFY(setupflag, "Call FindPointsGSLIB::Setup method first");
   auto *findptsData3 = (gslib::findpts_data_3 *)this->fdataD;
   auto *findptsData2 = (gslib::findpts_data_2 *)this->fdataD;
   int nve   = spacedim == 2 ? 4 : 8;
   int nel = NE_split_total;
   aabb.SetSize(spacedim*nve*nel);
   auto h_bb_ptr = DEV.bb.HostRead();
 
   if (spacedim == 3)
   {
      for (int e = 0; e < nel; e++)
      {
         Vector minn(spacedim), maxx(spacedim);
         if (dim == spacedim)
         {
            auto box = findptsData3->local.obb[e];
            for (int d = 0; d < spacedim; d++)
            {
               minn[d] = box.x[d].min;
               maxx[d] = box.x[d].max;
            }
         }
         else
         {
            const int n_el_ents = obb_check ?
                                  (3*spacedim + spacedim*spacedim) : (2*spacedim);
            const int min_off = obb_check ? spacedim : 0;
            const int max_off = obb_check ? 2*spacedim : spacedim;
            for (int d = 0; d < spacedim; d++)
            {
               minn[d] = h_bb_ptr[e*n_el_ents + min_off + d];
               maxx[d] = h_bb_ptr[e*n_el_ents + max_off + d];
            }
         }
         int c = 0;
         aabb(e*nve*spacedim + c++) = minn[0]; /* first vertex - x */
         aabb(e*nve*spacedim + c++) = minn[1]; /* y */
         aabb(e*nve*spacedim + c++) = minn[2]; /* z */
         aabb(e*nve*spacedim + c++) = maxx[0]; /* second vertex - x */
         aabb(e*nve*spacedim + c++) = minn[1]; /* . */
         aabb(e*nve*spacedim + c++) = minn[2]; /* . */
         aabb(e*nve*spacedim + c++) = maxx[0];
         aabb(e*nve*spacedim + c++) = maxx[1];
         aabb(e*nve*spacedim + c++) = minn[2];
         aabb(e*nve*spacedim + c++) = minn[0];
         aabb(e*nve*spacedim + c++) = maxx[1];
         aabb(e*nve*spacedim + c++) = minn[2];
         aabb(e*nve*spacedim + c++) = minn[0];
         aabb(e*nve*spacedim + c++) = minn[1];
         aabb(e*nve*spacedim + c++) = maxx[2];
         aabb(e*nve*spacedim + c++) = maxx[0];
         aabb(e*nve*spacedim + c++) = minn[1];
         aabb(e*nve*spacedim + c++) = maxx[2];
         aabb(e*nve*spacedim + c++) = maxx[0];
         aabb(e*nve*spacedim + c++) = maxx[1];
         aabb(e*nve*spacedim + c++) = maxx[2];
         aabb(e*nve*spacedim + c++) = minn[0];
         aabb(e*nve*spacedim + c++) = maxx[1];
         aabb(e*nve*spacedim + c++) = maxx[2];
      }
   }
   else // spacedim = 2
   {
      for (int e = 0; e < nel; e++)
      {
         Vector minn(spacedim), maxx(spacedim);
         if (dim == spacedim)
         {
            auto box = findptsData2->local.obb[e];
            for (int d = 0; d < spacedim; d++)
            {
               minn[d] = box.x[d].min;
               maxx[d] = box.x[d].max;
            }
         }
         else
         {
            const int n_el_ents = obb_check ? (3*spacedim + spacedim*spacedim) :
                                  (2*spacedim);
            const int min_off = obb_check ? spacedim : 0;
            const int max_off = obb_check ? 2*spacedim : spacedim;
            for (int d = 0; d < spacedim; d++)
            {
               minn[d] = h_bb_ptr[e*n_el_ents + min_off + d];
               maxx[d] = h_bb_ptr[e*n_el_ents + max_off + d];
            }
         }
         aabb(e*nve*spacedim + 0) = minn[0]; /* first vertex - x */
         aabb(e*nve*spacedim + 1) = minn[1]; /* y */
         aabb(e*nve*spacedim + 2) = maxx[0]; /* second vertex - x */
         aabb(e*nve*spacedim + 3) = minn[1]; /* . */
         aabb(e*nve*spacedim + 4) = maxx[0]; /* . */
         aabb(e*nve*spacedim + 5) = maxx[1];
         aabb(e*nve*spacedim + 6) = minn[0];
         aabb(e*nve*spacedim + 7) = maxx[1];
      }
   }
}
 
Mesh* FindPointsGSLIB::GetBoundingBoxMesh(int type)
{
   MFEM_VERIFY(setupflag, "Call FindPointsGSLIB::Setup method first");
   if (type != 0)
   {
      MFEM_VERIFY(obb_check || dim == spacedim,
                  "Oriented bounding boxes are not available when obb_check is false");
   }
   const unsigned int save_rank = 0;
   const unsigned int myid = gsl_comm->id;
   Vector bbvert;
   if (type == 0)
   {
      GetAxisAlignedBoundingBoxes(bbvert);
   }
   else
   {
      DenseTensor obbA;
      Vector obbC;
      GetOrientedBoundingBoxes(obbA, obbC, bbvert);
   }
   int nve   = spacedim == 2 ? 4 : 8;
   int nel = NE_split_total;
   int ne_glob = nel;
#ifdef MFEM_USE_MPI
   MPI_Allreduce(&nel, &ne_glob, 1, MPI_INT, MPI_SUM, gsl_comm->c);
#endif
 
 
   int nverts = nve*ne_glob;
   Mesh *meshbb = nullptr;
   if (myid == save_rank)
   {
      meshbb = new Mesh(spacedim, nverts, ne_glob, 0, spacedim);
   }
 
   int nsend = nel*nve*spacedim;
   MFEM_VERIFY(nsend == bbvert.Size(),
               "Inconsistent size of bounding box vertices");
   int nrecv = 0;
#ifdef MFEM_USE_MPI
   MPI_Status status;
#endif
   int vidx = 0;
   int eidx = 0;
   if (myid == save_rank)
   {
      for (int p = 0; (unsigned)p < gsl_comm->np; p++)
      {
         if (static_cast<unsigned int>(p) != save_rank)
         {
#ifdef MFEM_USE_MPI
            MPI_Recv(&nrecv, 1, MPI_INT, p, 444, gsl_comm->c, &status);
            bbvert.SetSize(nrecv);
            if (nrecv)
            {
               MPI_Recv(bbvert.GetData(), nrecv, MPI_DOUBLE, p, 445, gsl_comm->c, &status);
            }
#endif
         }
         else
         {
            nrecv = nsend;
         }
         int nel_recv = nrecv/(spacedim*nve);
         for (int e = 0; e < nel_recv; e++)
         {
            for (int j = 0; j < nve; j++)
            {
               Vector ver(bbvert.GetData() + e*nve*spacedim + j*spacedim, spacedim);
               meshbb->AddVertex(ver);
            }
 
            if (spacedim == 2)
            {
               const int inds[4] = {vidx++, vidx++, vidx++, vidx++};
               int attr = eidx+1;
               meshbb->AddQuad(inds, attr);
               eidx++;
            }
            else
            {
               const int inds[8] = {vidx++, vidx++, vidx++, vidx++,
                                    vidx++, vidx++, vidx++, vidx++
                                   };
               meshbb->AddHex(inds, (eidx++)+1);
            }
         }
      }
      if (spacedim == 2)
      {
         meshbb->FinalizeQuadMesh(1, 1, true);
      }
      else
      {
         meshbb->FinalizeHexMesh(1, 1, true);
      }
   }
   else
   {
#ifdef MFEM_USE_MPI
      MPI_Send(&nsend, 1, MPI_INT, save_rank, 444, gsl_comm->c);
      if (nsend)
      {
         MPI_Send(bbvert.GetData(), nsend, MPI_DOUBLE, save_rank, 445, gsl_comm->c);
      }
#endif
   }
#ifdef MFEM_USE_MPI
   MPI_Barrier(gsl_comm->c);
#endif
 
   return meshbb;
}
 
void FindPointsGSLIB::GetOrientedBoundingBoxes(DenseTensor &obbA, Vector &obbC,
                                               Vector &obbV) const
{
   MFEM_VERIFY(setupflag, "Call FindPointsGSLIB::Setup method first");
   MFEM_VERIFY(obb_check || dim == spacedim,
               "Oriented bounding boxes are not available when obb_check is false");
   auto *findptsData3 = (gslib::findpts_data_3 *)this->fdataD;
   auto *findptsData2 = (gslib::findpts_data_2 *)this->fdataD;
   int nve   = spacedim == 2 ? 4 : 8;
   int nel = NE_split_total;
 
   obbA.SetSize(spacedim, spacedim, nel);
   obbC.SetSize(spacedim*nel);
   obbV.SetSize(spacedim*nve*nel);
   if (spacedim == 3)
   {
      for (int e = 0; e < nel; e++)
      {
         double *Ad = obbA.GetData(e);
         if (dim == spacedim)
         {
            auto box = findptsData3->local.obb[e];
            for (int d = 0; d < spacedim; d++)
            {
               obbC(e*spacedim + d) = box.c0[d];
            }
            for (int i = 0; i < spacedim; i++)
            {
               for (int j = 0; j < spacedim; j++)
               {
                  Ad[i*spacedim + j] = box.A[i + j*spacedim]; // GSLIB uses row-major storage
               }
            }
         }
         else
         {
            int n_el_ents = 18;
            for (int d = 0; d < spacedim; d++)
            {
               obbC(e*spacedim + d) = DEV.bb(n_el_ents*e + d);
            }
            for (int i = 0; i < spacedim; i++)
            {
               for (int j = 0; j < spacedim; j++)
               {
                  Ad[i*spacedim + j] = DEV.bb[n_el_ents*e + 9 + j*spacedim+i];
               }
            }
         }
 
         DenseMatrix Amat = obbA(e);
         Amat.Invert();
         Vector center(obbC.GetData() + e*spacedim, spacedim);
 
         Vector v1(spacedim);
         Vector temp;
         v1(0) = -1.0; v1(1) = -1.0; v1(2) = -1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 0, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
         v1(0) = 1.0; v1(1) = -1.0; v1(2) = -1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 3, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
         v1(0) = 1.0; v1(1) = 1.0; v1(2) = -1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 6, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
         v1(0) = -1.0; v1(1) = 1.0; v1(2) = -1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 9, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
         v1(0) = -1.0; v1(1) = -1.0; v1(2) = 1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 12, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
         v1(0) = 1.0; v1(1) = -1.0; v1(2) = 1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 15, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
         v1(0) = 1.0; v1(1) = 1.0; v1(2) = 1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 18, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
         v1(0) = -1.0; v1(1) = 1.0; v1(2) = 1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 21, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
      }
   }
   else // spacedim = 2
   {
      for (int e = 0; e < nel; e++)
      {
         double *Ad = obbA.GetData(e);
         if (dim == spacedim)
         {
            auto box = findptsData2->local.obb[e];
            for (int d = 0; d < spacedim; d++)
            {
               obbC(e*spacedim + d) = box.c0[d];
            }
            for (int i = 0; i < spacedim; i++)
            {
               for (int j = 0; j < spacedim; j++)
               {
                  Ad[i*spacedim + j] = box.A[i + j*spacedim]; // GSLIB uses row-major storage
               }
            }
         }
         else
         {
            int n_el_ents = 10;
            for (int d = 0; d < spacedim; d++)
            {
               obbC(e*spacedim + d) = DEV.bb(n_el_ents*e + d);
            }
            for (int i = 0; i < spacedim; i++)
            {
               for (int j = 0; j < spacedim; j++)
               {
                  Ad[i*spacedim + j] = DEV.bb[n_el_ents*e + 6 + j*spacedim+i];
               }
            }
         }
 
 
         DenseMatrix Amat = obbA(e);
         Amat.Invert();
         Vector center(obbC.GetData() + e*spacedim, spacedim);
 
         Vector v1(spacedim);
         Vector temp;
         v1(0) = -1.0; v1(1) = -1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 0, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
         v1(0) = 1.0; v1(1) = -1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 2, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
         v1(0) = 1.0; v1(1) = 1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 4, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
         v1(0) = -1.0; v1(1) = 1.0;
         temp.SetDataAndSize(obbV.GetData() + e*nve*spacedim + 6, spacedim);
         Amat.Mult(v1, temp);
         temp += center;
      }
   }
}
 
void OversetFindPointsGSLIB::Setup(Mesh &m, const int meshid,
                                   GridFunction *gfmax,
                                   const double bbox_rel_size_inc,
                                   const double newt_tol,
                                   const int npt_max)
{
   MFEM_VERIFY(m.GetNodes() != nullptr, "Mesh nodes are required.");
   const int meshOrder = m.GetNodes()->FESpace()->GetMaxElementOrder();
   const int gfOrder = gfmax ? gfmax->FESpace()->GetMaxElementOrder() :
                       meshOrder;
   MFEM_VERIFY(meshOrder == gfOrder,
               "Mesh order must match gfmax order in OversetFindPointsGSLIB.");
 
   // FreeData if OversetFindPointsGSLIB::Setup has been called already
   if (setupflag) { FreeData(); }
   SetupCrystal();
 
   mesh = &m;
   dim  = mesh->Dimension();
   spacedim = dim;
   const FiniteElement *fe = mesh->GetNodalFESpace()->GetTypicalFE();
   unsigned dof1D = fe->GetOrder() + 1;
 
   SetupSplitMeshesAndIntegrationRules(meshOrder);
 
   GetNodalValues(mesh->GetNodes(), gsl_mesh);
 
   MFEM_ASSERT(meshid>=0, " The ID should be greater than or equal to 0.");
 
   const int pts_cnt = gsl_mesh.Size()/dim,
             NEtot = pts_cnt/(int)pow(dof1D, dim);
 
   distfint.SetSize(pts_cnt);
   if (!gfmax)
   {
      distfint = 0.0;
   }
   else
   {
      GetNodalValues(gfmax, distfint);
   }
   u_meshid = (unsigned int)meshid;
 
   if (dim == 2)
   {
      unsigned nr[2] = { dof1D, dof1D };
      unsigned mr[2] = { 2*dof1D, 2*dof1D };
      double * const elx[2] =
      {
         pts_cnt == 0 ? nullptr : &gsl_mesh(0),
         pts_cnt == 0 ? nullptr : &gsl_mesh(pts_cnt)
      };
      fdataD = findptsms_setup_2(gsl_comm, elx, nr, NEtot, mr,
                                 bbox_rel_size_inc, pts_cnt, pts_cnt,
                                 npt_max, newt_tol,
                                 &u_meshid, &distfint(0));
   }
   else
   {
      unsigned nr[3] = { dof1D, dof1D, dof1D };
      unsigned mr[3] = { 2*dof1D, 2*dof1D, 2*dof1D };
      double * const elx[3] =
      {
         pts_cnt == 0 ? nullptr : &gsl_mesh(0),
         pts_cnt == 0 ? nullptr : &gsl_mesh(pts_cnt),
         pts_cnt == 0 ? nullptr : &gsl_mesh(2*pts_cnt)
      };
      fdataD = findptsms_setup_3(gsl_comm, elx, nr, NEtot, mr,
                                 bbox_rel_size_inc, pts_cnt, pts_cnt,
                                 npt_max, newt_tol,
                                 &u_meshid, &distfint(0));
   }
   setupflag = true;
   overset   = true;
}
 
void OversetFindPointsGSLIB::FindPoints(const Vector &point_pos,
                                        const Array<unsigned int> &point_id,
                                        const int point_pos_ordering)
{
   MFEM_VERIFY(setupflag, "Use OversetFindPointsGSLIB::Setup before "
               "finding points.");
   MFEM_VERIFY(overset, "Please use OversetFindPoints for overlapping grids.");
   points_cnt = point_pos.Size() / dim;
   unsigned int match = 0; // Don't find points in the mesh if point_id=mesh_id
 
   gsl_code.SetSize(points_cnt);
   gsl_proc.SetSize(points_cnt);
   gsl_elem.SetSize(points_cnt);
   gsl_ref.SetSize(points_cnt * dim);
   gsl_dist.SetSize(points_cnt);
 
   auto xvFill = [&](const double *xv_base[], unsigned xv_stride[])
   {
      for (int d = 0; d < dim; d++)
      {
         if (point_pos_ordering == Ordering::byNODES)
         {
            xv_base[d] = point_pos.GetData() + d*points_cnt;
            xv_stride[d] = sizeof(double);
         }
         else
         {
            xv_base[d] = point_pos.GetData() + d;
            xv_stride[d] = dim*sizeof(double);
         }
      }
   };
   if (dim == 2)
   {
      auto *findptsData = (gslib::findpts_data_2 *)this->fdataD;
      const double *xv_base[2];
      unsigned xv_stride[2];
      xvFill(xv_base, xv_stride);
      findptsms_2(gsl_code.GetData(), sizeof(unsigned int),
                  gsl_proc.GetData(), sizeof(unsigned int),
                  gsl_elem.GetData(), sizeof(unsigned int),
                  gsl_ref.GetData(),  sizeof(double) * dim,
                  gsl_dist.GetData(), sizeof(double),
                  xv_base,            xv_stride,
                  point_id.GetData(), sizeof(unsigned int), &match,
                  points_cnt, findptsData);
   }
   else  // dim == 3
   {
      auto *findptsData = (gslib::findpts_data_3 *)this->fdataD;
      const double *xv_base[3];
      unsigned xv_stride[3];
      xvFill(xv_base, xv_stride);
      findptsms_3(gsl_code.GetData(), sizeof(unsigned int),
                  gsl_proc.GetData(), sizeof(unsigned int),
                  gsl_elem.GetData(), sizeof(unsigned int),
                  gsl_ref.GetData(),  sizeof(double) * dim,
                  gsl_dist.GetData(), sizeof(double),
                  xv_base,            xv_stride,
                  point_id.GetData(), sizeof(unsigned int), &match,
                  points_cnt, findptsData);
   }
 
   // Set the element number and reference position to 0 for points not found
   for (int i = 0; i < points_cnt; i++)
   {
      if (gsl_code[i] == 2 ||
          (gsl_code[i] == 1 && gsl_dist(i) > bdr_tol))
      {
         gsl_elem[i] = 0;
         for (int d = 0; d < dim; d++) { gsl_ref(i*dim + d) = -1.; }
         gsl_code[i] = 2;
      }
   }
 
   // Map element number for simplices, and ref_pos from [-1,1] to [0,1] for both
   // simplices and quads.
   MapRefPosAndElemIndices();
}
 
void OversetFindPointsGSLIB::Interpolate(const Vector &point_pos,
                                         const Array<unsigned int> &point_id,
                                         const GridFunction &field_in,
                                         Vector &field_out,
                                         const int point_pos_ordering)
{
   FindPoints(point_pos, point_id, point_pos_ordering);
   Interpolate(field_in, field_out);
}
 
GSOPGSLIB::GSOPGSLIB(Array<long long> &ids)
{
   gsl_comm = new gslib::comm;
   cr       = new gslib::crystal;
#ifdef MFEM_USE_MPI
   if (!Mpi::IsInitialized()) { Mpi::Init(); }
   MPI_Comm comm = MPI_COMM_WORLD;
   comm_init(gsl_comm, comm);
#else
   comm_init(gsl_comm, 0);
#endif
   crystal_init(cr, gsl_comm);
   UpdateIdentifiers(ids);
}
 
#ifdef MFEM_USE_MPI
GSOPGSLIB::GSOPGSLIB(MPI_Comm comm_, Array<long long> &ids)
   : cr(NULL), gsl_comm(NULL)
{
   gsl_comm = new gslib::comm;
   cr      = new gslib::crystal;
   comm_init(gsl_comm, comm_);
   crystal_init(cr, gsl_comm);
   UpdateIdentifiers(ids);
}
#endif
 
GSOPGSLIB::~GSOPGSLIB()
{
#ifdef MFEM_USE_MPI
   if (!Mpi::IsFinalized())  // currently segfaults inside gslib otherwise
#endif
   {
      crystal_free(cr);
      comm_free(gsl_comm);
      delete gsl_comm;
      delete cr;
      gslib_gs_free(gsl_data);
   }
}
 
void GSOPGSLIB::UpdateIdentifiers(const Array<long long> &ids)
{
   long long minval = ids.Min();
#ifdef MFEM_USE_MPI
   MPI_Allreduce(MPI_IN_PLACE, &minval, 1, MPI_LONG_LONG_INT,
                 MPI_MIN, gsl_comm->c);
#endif
   MFEM_VERIFY(minval >= 0, "Unique identifier cannot be negative.");
   if (gsl_data != nullptr) { gslib_gs_free(gsl_data); }
   num_ids = ids.Size();
   gsl_data = gslib_gs_setup(ids.GetData(),
                             ids.Size(),
                             gsl_comm, 0,
                             gslib::gs_crystal_router, 0);
}
 
void GSOPGSLIB::GS(Vector &senddata, GSOp op)
{
   MFEM_VERIFY(senddata.Size() == num_ids,
               "Incompatible setup and GOP operation.");
   if (op == GSOp::ADD)
   {
      gslib_gs(senddata.GetData(),gslib::gs_double,gslib::gs_add,0,gsl_data,0);
   }
   else if (op == GSOp::MUL)
   {
      gslib_gs(senddata.GetData(),gslib::gs_double,gslib::gs_mul,0,gsl_data,0);
   }
   else if (op == GSOp::MAX)
   {
      gslib_gs(senddata.GetData(),gslib::gs_double,gslib::gs_max,0,gsl_data,0);
   }
   else if (op == GSOp::MIN)
   {
      gslib_gs(senddata.GetData(),gslib::gs_double,gslib::gs_min,0,gsl_data,0);
   }
   else
   {
      MFEM_ABORT("Invalid GSOp operation.");
   }
}
 
#if defined(MFEM_USE_MPI)
void GlobalBBoxTensorGridMap::SetupCrystal(const MPI_Comm &comm_)
{
#if GSLIB_RELEASE_VERSION < 10009
   MFEM_ABORT("GSLIB version 1.0.9 or higher is required.");
#endif
   gsl_comm = new gslib::comm;
   cr      = new gslib::crystal;
   comm_init(gsl_comm, comm_);
   gslib::crystal_init(cr, gsl_comm);
}
 
GlobalBBoxTensorGridMap::GlobalBBoxTensorGridMap(ParMesh &pmesh, int nx)
{
   GridFunction *nodes = pmesh.GetNodes();
   const int nel = pmesh.GetNE();
   sdim = pmesh.SpaceDimension();
   Vector elmin(nel*sdim), elmax(nel*sdim);
   elmin = std::numeric_limits<real_t>::max();
   elmax = -std::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++)
      {
         pmesh.GetElementVertices(e, verts);
         Vector center(sdim);
         for (int v = 0; v < verts.Size(); v++)
         {
            coord = pmesh.GetVertex(verts[v]);
            for (int d = 0; d < sdim; d++)
            {
               elmin(d*nel + e) = std::min(elmin(d*nel + e), coord[d]);
               elmax(d*nel + e) = std::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(pmesh.GetComm(), elmin, elmax, nel, nx_arr);
}
 
GlobalBBoxTensorGridMap::GlobalBBoxTensorGridMap(const MPI_Comm &comm,
                                                 Vector &elmin, Vector &elmax,
                                                 int nel, int sdim_,
                                                 int n,
                                                 bool by_max_size)
{
   sdim = sdim_;
   MFEM_VERIFY(0 < sdim && sdim <= 3,
               "GlobalBBoxTensorGridMap 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);
   if (!by_max_size)
   {
      nx_arr = n;
   }
   else
   {
      long long int nx = n;
      MPI_Allreduce(MPI_IN_PLACE, &nx, 1, MPI_LONG_LONG, MPI_SUM, comm);
      nx = ceil(pow((double)nx,1./sdim));
      nx_arr = nx;
   }
   Setup(comm, elmin, elmax, nel, nx_arr);
}
 
GlobalBBoxTensorGridMap::GlobalBBoxTensorGridMap(const MPI_Comm &comm,
                                                 Vector &elmin, Vector &elmax,
                                                 int nel, int sdim_,
                                                 Array<int> &nx)
{
   sdim = sdim_;
   Setup(comm, elmin, elmax, nel, nx);
}
 
void GlobalBBoxTensorGridMap::Setup(const MPI_Comm &comm,
                                    Vector &elmin, Vector &elmax,
                                    int nel, Array<int> &nx)
{
   SetupCrystal(comm);
   MFEM_VERIFY(0 < sdim && sdim <= 3,
               "GlobalBBoxTensorGridMap only supports spatial dimensions 1, 2, and 3.");
   MFEM_VERIFY(nx.Size() == sdim,
               "GlobalBBoxTensorGridMap 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.");
   }
   gmap_bnd_min.SetSize(sdim);
   gmap_bnd_max.SetSize(sdim);
   gmap_fac.SetSize(sdim);
   gmap_n.SetSize(sdim);
   gmap_n = nx;
 
   long long int global_nel = nel;
   MPI_Allreduce(MPI_IN_PLACE, &global_nel, 1, MPI_LONG_LONG, MPI_SUM, comm);
 
   MPI_Comm_size(comm, &num_procs);
   int gmap_nd = gmap_n[0];
   for (int d = 1; d < sdim; d++)
   {
      gmap_nd *= gmap_n[d];
   }
 
   if (global_nel == 0)
   {
      gmap_n = 1;
      gmap_nd = 1;
 
      // Mark the global bounding box as empty so all point queries return no
      // candidate ranks without relying on sentinel infinities.
      gmap_bnd_min = 1.0;
      gmap_bnd_max = 0.0;
      gmap_fac = 0.0;
      n_local_cells = (gmap_nd - 1) / num_procs + 1;
      ggrid_map.SetSize(n_local_cells + 1);
      ggrid_map = n_local_cells + 1;
      return;
   }
 
   for (int d = 0; d < nx.Size(); d++)
   {
      MFEM_VERIFY(nx[d] > 0,
                  "GlobalBBoxTensorGridMap requires positive number of divisions in each dimension.");
   }
   gmap_bnd_min = std::numeric_limits<real_t>::max();
   gmap_bnd_max = -std::numeric_limits<real_t>::max();
   if (nel > 0)
   {
      for (int d = 0; d < sdim; d++)
      {
         Vector elmind(elmin.GetData() + d*nel, nel);
         Vector elmaxd(elmax.GetData() + d*nel, nel);
         gmap_bnd_min[d] = elmind.Min();
         gmap_bnd_max[d] = elmaxd.Max();
      }
   }
 
   Vector gmap_bnd_min_loc = gmap_bnd_min;
   Vector gmap_bnd_max_loc = gmap_bnd_max;
 
   MPI_Allreduce(MPI_IN_PLACE, gmap_bnd_min.GetData(), sdim,
                 MFEM_MPI_REAL_T, MPI_MIN, comm);
   MPI_Allreduce(MPI_IN_PLACE, gmap_bnd_max.GetData(), sdim,
                 MFEM_MPI_REAL_T, MPI_MAX, comm);
 
   BBoxTensorGridMap::SetGridFac(gmap_fac, gmap_n, gmap_bnd_min, gmap_bnd_max);
 
   // Grid cell ranges for each element in each dimension
   Array<int> elmin_h, elmax_h;
   int store_size = BBoxTensorGridMap::GetGridCountAndRange(gmap_n, gmap_fac,
                                                            gmap_bnd_min,
                                                            gmap_bnd_max,
                                                            elmin, elmax,
                                                            elmin_h, elmax_h);
 
   Array<int> loc_idx_min(sdim), loc_idx_max(sdim), lh_n(sdim);
   int loc_idx_tot = 1;
 
   if (nel > 0)
   {
      for (int d = 0; d < sdim; d++)
      {
         BBoxTensorGridMap::GetGridRange(d, gmap_n, gmap_fac,
                                         gmap_bnd_min,
                                         gmap_bnd_min_loc[d],
                                         gmap_bnd_max_loc[d],
                                         loc_idx_min[d], loc_idx_max[d]);
         lh_n[d] = loc_idx_max[d] - loc_idx_min[d];
         loc_idx_tot *= lh_n[d];
      }
   }
   else
   {
      loc_idx_min = 0;
      loc_idx_max = 1;
      lh_n = 1;
   }
 
   struct hashInfo_s
   {
      unsigned int index, proc;
   };
   struct gslib::array hashInfo_pt;
   array_init(struct hashInfo_s, &hashInfo_pt, store_size);
   hashInfo_pt.n=store_size;
 
   struct hashInfo_s *pt = (struct hashInfo_s *)hashInfo_pt.ptr;
   Array<int> marker(loc_idx_tot);
   marker = 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) * gmap_n[0] * gmap_n[1];
         int koff_loc = sdim < 3 ? 0 :
                        (elmin_h[2*nel+e]+k - loc_idx_min[2])*lh_n[0] * lh_n[1];
         for (int j = 0; j < jlim; j++)
         {
            int joff = sdim < 2 ? 0 : (elmin_h[1*nel + e] + j) * gmap_n[0];
            int joff_loc = sdim < 2 ? 0 :
                           (elmin_h[1*nel+e]+j - loc_idx_min[1])*lh_n[0];
            for (int i = 0; i < ilim; i++)
            {
               int ioff = elmin_h[0*nel + e] + i;
               int ioff_loc = elmin_h[0*nel+e]+i - loc_idx_min[0];
               int idx = ioff + joff + koff;
               int idx_loc = ioff_loc + joff_loc + koff_loc;
               if (marker[idx_loc] == 1) { continue; }
               pt->proc = idx % num_procs;
               pt->index = idx / num_procs;
               marker[idx_loc] = 1;
               ++pt;
            }
         }
      }
   }
   MPI_Barrier(comm);
   int npts = marker.Sum();
   hashInfo_pt.n = npts;
 
   // transfer info to other ranks and sort by index
#if GSLIB_RELEASE_VERSION >= 10009
   sarray_transfer(struct hashInfo_s, &hashInfo_pt, proc, 1, cr);
   sarray_sort(struct hashInfo_s, hashInfo_pt.ptr, hashInfo_pt.n,
               index, 0, &(cr->data));
#endif
 
   int nrecv = hashInfo_pt.n;
 
   n_local_cells = (gmap_nd-1)/num_procs+1;
   ggrid_map.SetSize(n_local_cells + 1 + nrecv);
   ggrid_map = 0;
   Array<int> hash_el_count(n_local_cells);
   hash_el_count = 0;
 
   pt = (struct hashInfo_s *)hashInfo_pt.ptr;
   for (int i = 0; i < nrecv; i++)
   {
      int idx = pt->index;
      hash_el_count[idx]++;
      ++pt;
   }
 
   ggrid_map[0] = n_local_cells + 1;
   for (int e = 0; e < n_local_cells; e++)
   {
      ggrid_map[e + 1] = ggrid_map[e] + hash_el_count[e];
   }
 
   pt = (struct hashInfo_s *)hashInfo_pt.ptr;
   for (int i = 0; i < nrecv; i++)
   {
      int idx = pt->index;
      int proc = pt->proc;
      ggrid_map[ggrid_map[idx+1]-hash_el_count[idx]]=proc;
      hash_el_count[idx]--;
      ++pt;
   }
 
   array_free(&hashInfo_pt);
   MPI_Barrier(comm);
}
 
int GlobalBBoxTensorGridMap::GetGlobalGridCellFromPoint(Vector &xyz) const
{
   MFEM_ASSERT(xyz.Size() == sdim,
               "Point must have the same dimension as the hash.");
   int sum = 0;
   for (int d = sdim-1; d >= 0; --d)
   {
      if (xyz(d) < gmap_bnd_min(d) || xyz(d) > gmap_bnd_max(d))
      {
         return -1; // Point is outside the bounds of the hash
      }
      sum *= gmap_n[d];
      int i = (int)floor((xyz(d) - gmap_bnd_min(d)) * gmap_fac[d]);
      sum += i < 0 ? 0 : (gmap_n[d] - 1 < i ? gmap_n[d] - 1 : i);
   }
   return sum;
}
 
void GlobalBBoxTensorGridMap::GlobalGridCellToProcAndLocalIndex(int i,
                                                                int &proc,
                                                                int &idx) const
{
   proc = i % num_procs;
   idx = i / num_procs;
}
 
void GlobalBBoxTensorGridMap::GetProcAndLocalIndexFromPoint(Vector &xyz,
                                                            int &proc, int &idx) const
{
   int cell = GetGlobalGridCellFromPoint(xyz);
   if (cell < 0)
   {
      proc = -1;
      idx = -1;
      return; // Point is outside the bounds of the grid
   }
   GlobalGridCellToProcAndLocalIndex(cell, proc, idx);
}
 
void GlobalBBoxTensorGridMap::MapPointsToProcs(Vector &xyz, int ordering,
                                               std::map<int,
                                               std::vector<int>> &pt_to_procs) const
{
   MFEM_ASSERT(xyz.Size() % sdim == 0,
               "Point array size must be a multiple of the grid dimension.");
   int npts = xyz.Size() / sdim;
   pt_to_procs.clear();
   for (int i = 0; i < npts; i++)
   {
      pt_to_procs[i]; // ensure key exists for every point
   }
 
   // struct to hold point info:
   // info: holds local grid cell index when first sent.
   //       holds proc when returned back.
   // proc: holds the proc that the grid cell info is on.
   // loc_index: local index of the point in the input vector.
   struct ptInfo_s
   {
      unsigned int info, proc, loc_index;
   };
   struct gslib::array ptInfo_pt;
   array_init(struct ptInfo_s, &ptInfo_pt, npts);
   struct ptInfo_s *pt = (struct ptInfo_s *)ptInfo_pt.ptr;
   struct ptInfo_s *pt_begin = pt;
 
   for (int i = 0; i < npts; i++)
   {
      Vector pt_xyz(sdim);
      for (int d = 0; d < sdim; d++)
      {
         pt_xyz(d) = ordering == 0 ? xyz(d*npts + i) : xyz(i*sdim + d);
      }
 
      int cell = GetGlobalGridCellFromPoint(pt_xyz);
      if (cell < 0)
      {
         continue; // Point is outside the bounds of the hash
      }
      int proc, idx;
      GlobalGridCellToProcAndLocalIndex(cell, proc, idx);
      pt->info = idx; // local grid cell index on proc
      pt->proc = proc; // proc that the grid cell info is on
      pt->loc_index = i; // local index of the point in the input vector
      ++pt;
   }
   ptInfo_pt.n = pt - pt_begin;
   MPI_Barrier(gsl_comm->c);
 
   // transfer info to ranks that hold each hash cell's info
   sarray_transfer(struct ptInfo_s, &ptInfo_pt, proc, 1, cr);
 
   int nrecv = ptInfo_pt.n;
   pt = (struct ptInfo_s *)ptInfo_pt.ptr;
   int ncount = 0;
   for (int i = 0; i < nrecv; i++)
   {
      int idx = pt->info;
      int loc_count = ggrid_map[idx+1]-ggrid_map[idx];
      ncount += loc_count;
      ++pt;
   }
 
   struct gslib::array sendptInfo_pt;
   array_init(struct ptInfo_s, &sendptInfo_pt, ncount);
   sendptInfo_pt.n=ncount;
   struct ptInfo_s *spt = (struct ptInfo_s *)sendptInfo_pt.ptr;
   pt = (struct ptInfo_s *)ptInfo_pt.ptr;
 
   for (int i = 0; i < nrecv; i++)
   {
      Array<int> procs = MapCellToProcs(pt->info); // map for local index
      for (int k = 0; k < procs.Size(); k++)
      {
         spt->info = procs[k]; // proc that has contribution to this cell
         spt->proc = pt->proc; // proc that sent this point
         spt->loc_index = pt->loc_index; // local index of the point originally
         ++spt;
      }
      ++pt;
   }
 
   array_free(&ptInfo_pt);
   // transfer information back to ranks where the query originated from.
   sarray_transfer(struct ptInfo_s, &sendptInfo_pt, proc, 1, cr);
 
   nrecv = sendptInfo_pt.n;
   spt = (struct ptInfo_s *)sendptInfo_pt.ptr;
 
   for (int i =0; i < nrecv; i++)
   {
      int pt_idx = spt->loc_index; // local index of the point
      int proc = spt->info; // proc that could be holding this point
      pt_to_procs[pt_idx].push_back(proc);
      ++spt;
   }
 
   array_free(&sendptInfo_pt);
   MPI_Barrier(gsl_comm->c);
}
 
Array<int> GlobalBBoxTensorGridMap::MapCellToProcs(int l_idx) const
{
   MFEM_ASSERT(l_idx >= 0 && l_idx < n_local_cells,
               "Access element " << l_idx << " of local hash with cells = "
               << n_local_cells);
   int start = ggrid_map[l_idx];
   int end = ggrid_map[l_idx + 1];
   Array<int> elements(end - start);
   for (int j = start; j < end; j++)
   {
      elements[j - start] = ggrid_map[j];
   }
   return elements;
}
 
GlobalBBoxTensorGridMap::~GlobalBBoxTensorGridMap()
{
   if (!Mpi::IsFinalized())  // currently segfaults inside gslib otherwise
   {
      crystal_free(cr);
      comm_free(gsl_comm);
      delete gsl_comm;
      delete cr;
   }
}
 
#endif // defined(MFEM_USE_MPI) && defined(MFEM_USE_GSLIB)
 
} // namespace mfem
#undef CODE_INTERNAL
#undef CODE_BORDER
#undef CODE_NOT_FOUND
 
#endif // MFEM_USE_GSLIB