interpolate_local_3.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 "../../general/forall.hpp"
#include "gslib_kernel_helpers.hpp"
 
#ifdef MFEM_USE_GSLIB
 
#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-function"
#endif
#include "gslib.h"
#ifndef GSLIB_RELEASE_VERSION //gslib v1.0.7
#define GSLIB_RELEASE_VERSION 10007
#endif
#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
namespace mfem
{
#if GSLIB_RELEASE_VERSION >= 10009
#define CODE_INTERNAL 0
#define CODE_BORDER 1
#define CODE_NOT_FOUND 2
 
using gslib::lagrange_eval;
 
template<int T_D1D = 0>
static void InterpolateLocal3DKernel(const double *const gf_in,
                                     int *const el,
                                     double *const r,
                                     double *const int_out,
                                     const int npt,
                                     const int ncomp,
                                     double *gll1D,
                                     double *lagcoeff,
                                     const int pN = 0)
{
   const int Nfields = ncomp;
   const int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
   const int D1D = T_D1D ? T_D1D : pN;
   const int p_Np = D1D*D1D*D1D;
   MFEM_VERIFY(MD1 <= DofQuadLimits::MAX_D1D,
               "Increase Max allowable polynomial order.");
   MFEM_VERIFY(D1D != 0, "Polynomial order not specified.");
#define MAXC(a, b) (((a) > (b)) ? (a) : (b))
   const int nThreadsy = MAXC(D1D, 3);
   mfem::forall_2D(npt, D1D, nThreadsy, [=] MFEM_HOST_DEVICE (int i)
   {
      MFEM_SHARED double wtr[3*MD1];
      MFEM_SHARED double sums[MD1*MD1];
 
      // Evaluate basis functions at the reference space coordinates
      MFEM_FOREACH_THREAD(j,x,D1D)
      {
         MFEM_FOREACH_THREAD(k,y,3)
         {
            lagrange_eval(wtr + k*D1D, r[3*i+k], j, D1D, gll1D, lagcoeff);
         }
      }
      MFEM_SYNC_THREAD;
 
      for (int fld = 0; fld < Nfields; ++fld)
      {
         // If using GetNodalValues, ordering is NDOFS x NEL x VDIM and the
         // offset would be `el[i] * p_Np + fld * gf_offset`.
         // R->Mult produces element vectors in NDOFS x VDIM x NEL layout.
         const int elemOffset = el[i] * p_Np * Nfields + fld * p_Np;
         MFEM_FOREACH_THREAD(j,x,D1D)
         {
            MFEM_FOREACH_THREAD(k,y,D1D)
            {
               sums[j + k*D1D] = 0.0;
               for (int l = 0; l < D1D; ++l)
               {
                  sums[j + k*D1D] += gf_in[elemOffset + j + k*D1D + l*D1D*D1D] *
                                     wtr[2*D1D+l];
               }
               sums[j+k*D1D] *= wtr[D1D+k]*wtr[j];
            }
         }
         MFEM_SYNC_THREAD;
 
         MFEM_FOREACH_THREAD(j,x,1)
         {
            MFEM_FOREACH_THREAD(k,y,1)
            {
               double sumv = 0.0;
               for (int jj = 0; jj < D1D*D1D; ++jj)
               {
                  sumv += sums[jj];
               }
               int_out[i + fld * npt] = sumv;
            }
         }
         MFEM_SYNC_THREAD;
      }
   });
}
 
void FindPointsGSLIB::InterpolateLocal3(const Vector &field_in,
                                        Array<int> &gsl_elem_dev_l,
                                        Vector &gsl_ref_l,
                                        Vector &field_out,
                                        int npt, int ncomp,
                                        int dof1Dsol)
{
   if (npt == 0) { return; }
   bool use_dev = field_in.UseDevice();
   auto pfin = field_in.Read(use_dev);
   auto pgsle = gsl_elem_dev_l.ReadWrite(use_dev);
   auto pgslr = gsl_ref_l.ReadWrite(use_dev);
   auto pfout = field_out.Write(use_dev);
   auto pgll = DEV.gll1d_sol.ReadWrite(use_dev);
   auto plcf = DEV.lagcoeff_sol.ReadWrite(use_dev);
   switch (dof1Dsol)
   {
      case 2:
         InterpolateLocal3DKernel<2>(pfin, pgsle, pgslr, pfout,
                                     npt, ncomp, pgll, plcf);
         break;
      case 3:
         InterpolateLocal3DKernel<3>(pfin, pgsle, pgslr, pfout,
                                     npt, ncomp, pgll, plcf);
         break;
      case 4:
         InterpolateLocal3DKernel<4>(pfin, pgsle, pgslr, pfout,
                                     npt, ncomp, pgll, plcf);
         break;
      case 5:
         InterpolateLocal3DKernel<5>(pfin, pgsle, pgslr, pfout,
                                     npt, ncomp, pgll, plcf);
         break;
      default:
         InterpolateLocal3DKernel(pfin, pgsle, pgslr, pfout,
                                  npt, ncomp, pgll, plcf, dof1Dsol);
         break;
   }
}
 
 
#undef MAXC
#undef CODE_INTERNAL
#undef CODE_BORDER
#undef CODE_NOT_FOUND
#else
void FindPointsGSLIB::InterpolateLocal3(const Vector &field_in,
                                        Array<int> &gsl_elem_dev_l,
                                        Vector &gsl_ref_l,
                                        Vector &field_out,
                                        int npt, int ncomp,
                                        int dof1Dsol) {};
#endif
} // namespace mfem
 
#endif //ifdef MFEM_USE_GSLIB