Formatting.

CHANGELOG

@@ -56,6 +56,9 @@ Discretization improvements
 
 - Added support for p-refined meshes in FindPointsGSLIB.
 
+- Added radial basis functions and reproducing kernels as alternatives to 
+  polynomial functions. These are accessed through the LocalKernelFECollection.
+
 Linear and nonlinear solvers
 ----------------------------
 - Updated interface to MUMPS direct solver to support multiple right-hand

examples/ex14.cpp

@@ -18,6 +18,8 @@
 //               ex14 -m ../data/amr-quad.mesh -r 3
 //               ex14 -m ../data/amr-hex.mesh
 //               ex14 -m ../data/fichera-amr.mesh
+//               ex14 -m ../data/inline-quad.mesh -r 2 -rk
+//               ex14 -m ../data/star.mesh -r 0 -o 3 -rk
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               discontinuous Galerkin (DG) finite element discretization of
@@ -43,10 +45,11 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/star.mesh";
    int ref_levels = -1;
    int order = 1;
+   bool rk = false;
    double sigma = -1.0;
    double kappa = -1.0;
    double eta = 0.0;
-   bool visualization = 1;
+   bool visualization = true;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -54,7 +57,9 @@ int main(int argc, char *argv[])
    args.AddOption(&ref_levels, "-r", "--refine",
                   "Number of times to refine the mesh uniformly, -1 for auto.");
    args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree) >= 0.");
+                  "Finite element order (polynomial degree) >= 0 OR reproducing kernel order.");
+   args.AddOption(&rk, "-rk", "--rk", "-no-rk", "--no-rk",
+                  "Use reproducing kernel functions");
    args.AddOption(&sigma, "-s", "--sigma",
                   "One of the three DG penalty parameters, typically +1/-1."
                   " See the documentation of class DGDiffusionIntegrator.");
@@ -75,6 +80,10 @@ int main(int argc, char *argv[])
    {
       kappa = (order+1)*(order+1);
    }
+   if (rk && sigma < 0.0)
+   {
+      sigma = 1.0;
+   }
    args.PrintOptions(cout);
 
    // 2. Read the mesh from the given mesh file. We can handle triangular,
@@ -104,7 +113,15 @@ int main(int argc, char *argv[])
 
    // 4. Define a finite element space on the mesh. Here we use discontinuous
    //    finite elements of the specified order >= 0.
-   FiniteElementCollection *fec = new DG_FECollection(order, dim);
+   FiniteElementCollection *fec;
+   if (rk)
+   {
+      fec = new LocalKernelFECollection(dim, 4, 6, order, 1.01 + order, 0.0);
+   }
+   else
+   {
+      fec = new DG_FECollection(order, dim);
+   }
    FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
    cout << "Number of unknowns: " << fespace->GetVSize() << endl;
 

examples/ex14p.cpp

@@ -17,6 +17,8 @@
 //               mpirun -np 4 ex14p -m ../data/inline-segment.mesh -rs 5
 //               mpirun -np 4 ex14p -m ../data/amr-quad.mesh -rs 3
 //               mpirun -np 4 ex14p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex14p -m ../data/star.mesh -rs 1 -o 2 -rk
+//               mpirun -np 4 ex14p -m ../data/amr-quad.mesh -rs 2 -rk
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               discontinuous Galerkin (DG) finite element discretization of
@@ -81,6 +83,7 @@ int main(int argc, char *argv[])
    int ser_ref_levels = -1;
    int par_ref_levels = 2;
    int order = 1;
+   bool rk = false;
    double sigma = -1.0;
    double kappa = -1.0;
    double eta = 0.0;
@@ -96,6 +99,8 @@ int main(int argc, char *argv[])
                   "Number of times to refine the mesh uniformly in parallel.");
    args.AddOption(&order, "-o", "--order",
                   "Finite element order (polynomial degree) >= 0.");
+   args.AddOption(&rk, "-rk", "--rk", "-no-rk", "--no-rk",
+                  "Use reproducing kernel functions");
    args.AddOption(&sigma, "-s", "--sigma",
                   "One of the three DG penalty parameters, typically +1/-1."
                   " See the documentation of class DGDiffusionIntegrator.");
@@ -119,6 +124,10 @@ int main(int argc, char *argv[])
    {
       kappa = (order+1)*(order+1);
    }
+   if (rk && sigma < 0.0)
+   {
+      sigma = 1.0;
+   }
    if (myid == 0)
    {
       args.PrintOptions(cout);
@@ -163,7 +172,15 @@ int main(int argc, char *argv[])
 
    // 6. Define a parallel finite element space on the parallel mesh. Here we
    //    use discontinuous finite elements of the specified order >= 0.
-   FiniteElementCollection *fec = new DG_FECollection(order, dim);
+   FiniteElementCollection *fec;
+   if (rk)
+   {
+      fec = new LocalKernelFECollection(dim, 4, 6, order, 1.01 + order, 0.0);
+   }
+   else
+   {
+      fec = new DG_FECollection(order, dim);
+   }
    ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
    HYPRE_BigInt size = fespace->GlobalTrueVSize();
    if (myid == 0)

examples/ex18.cpp

@@ -9,6 +9,8 @@
 //       ex18 -p 1 -r 0 -o 5 -s 6
 //       ex18 -p 2 -r 1 -o 1 -s 3
 //       ex18 -p 2 -r 0 -o 3 -s 3
+//       ex18 -p 1 -r 1 -o 3 -s 4 -rk
+//       ex18 -p 2 -r 0 -o 3 -s 3 -rk
 //
 // Description:  This example code solves the compressible Euler system of
 //               equations, a model nonlinear hyperbolic PDE, with a
@@ -65,6 +67,7 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/periodic-square.mesh";
    int ref_levels = 1;
    int order = 3;
+   bool rk = false;
    int ode_solver_type = 4;
    double t_final = 2.0;
    double dt = -0.01;
@@ -84,6 +87,8 @@ int main(int argc, char *argv[])
                   "Number of times to refine the mesh uniformly.");
    args.AddOption(&order, "-o", "--order",
                   "Order (degree) of the finite elements.");
+   args.AddOption(&rk, "-rk", "--rk", "-no-rk", "--no-rk",
+                  "Use reproducing kernel functions");
    args.AddOption(&ode_solver_type, "-s", "--ode-solver",
                   "ODE solver: 1 - Forward Euler,\n\t"
                   "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
@@ -139,13 +144,21 @@ int main(int argc, char *argv[])
 
    // 5. Define the discontinuous DG finite element space of the given
    //    polynomial order on the refined mesh.
-   DG_FECollection fec(order, dim);
+   FiniteElementCollection *fec;
+   if (rk)
+   {
+      fec = new LocalKernelFECollection(dim, 5, 7, order, 2.01 + order, 1.0);
+   }
+   else
+   {
+      fec = new DG_FECollection(order, dim);
+   }
    // Finite element space for a scalar (thermodynamic quantity)
-   FiniteElementSpace fes(&mesh, &fec);
+   FiniteElementSpace fes(&mesh, fec);
    // Finite element space for a mesh-dim vector quantity (momentum)
-   FiniteElementSpace dfes(&mesh, &fec, dim, Ordering::byNODES);
+   FiniteElementSpace dfes(&mesh, fec, dim, Ordering::byNODES);
    // Finite element space for all variables together (total thermodynamic state)
-   FiniteElementSpace vfes(&mesh, &fec, num_equation, Ordering::byNODES);
+   FiniteElementSpace vfes(&mesh, fec, num_equation, Ordering::byNODES);
 
    // This example depends on this ordering of the space.
    MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
@@ -304,6 +317,7 @@ int main(int argc, char *argv[])
 
    // Free the used memory.
    delete ode_solver;
+   delete fec;
 
    return 0;
 }

examples/ex18p.cpp

@@ -9,6 +9,8 @@
 //       mpirun -np 4 ex18p -p 1 -rs 1 -rp 1 -o 5 -s 6
 //       mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 1 -s 3
 //       mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 3 -s 3
+//       mpirun -np 4 ex18p -p 1 -rs 1 -rp 0 -o 3 -s 4 -rk
+//       mpirun -np 4 ex18p -p 2 -rs 1 -rp 0 -o 3 -s 3 -rk
 //
 // Description:  This example code solves the compressible Euler system of
 //               equations, a model nonlinear hyperbolic PDE, with a
@@ -70,6 +72,7 @@ int main(int argc, char *argv[])
    int ser_ref_levels = 0;
    int par_ref_levels = 1;
    int order = 3;
+   bool rk = false;
    int ode_solver_type = 4;
    double t_final = 2.0;
    double dt = -0.01;
@@ -93,6 +96,8 @@ int main(int argc, char *argv[])
                   " partitioning.");
    args.AddOption(&order, "-o", "--order",
                   "Order (degree) of the finite elements.");
+   args.AddOption(&rk, "-rk", "--rk", "-no-rk", "--no-rk",
+                  "Use reproducing kernel functions");
    args.AddOption(&ode_solver_type, "-s", "--ode-solver",
                   "ODE solver: 1 - Forward Euler,\n\t"
                   "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
@@ -161,13 +166,21 @@ int main(int argc, char *argv[])
 
    // 7. Define the discontinuous DG finite element space of the given
    //    polynomial order on the refined mesh.
-   DG_FECollection fec(order, dim);
+   FiniteElementCollection *fec;
+   if (rk)
+   {
+      fec = new LocalKernelFECollection(dim, 5, 7, order, 2.01 + order, 1.0);
+   }
+   else
+   {
+      fec = new DG_FECollection(order, dim);
+   }
    // Finite element space for a scalar (thermodynamic quantity)
-   ParFiniteElementSpace fes(&pmesh, &fec);
+   ParFiniteElementSpace fes(&pmesh, fec);
    // Finite element space for a mesh-dim vector quantity (momentum)
-   ParFiniteElementSpace dfes(&pmesh, &fec, dim, Ordering::byNODES);
+   ParFiniteElementSpace dfes(&pmesh, fec, dim, Ordering::byNODES);
    // Finite element space for all variables together (total thermodynamic state)
-   ParFiniteElementSpace vfes(&pmesh, &fec, num_equation, Ordering::byNODES);
+   ParFiniteElementSpace vfes(&pmesh, fec, num_equation, Ordering::byNODES);
 
    // This example depends on this ordering of the space.
    MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
@@ -375,6 +388,7 @@ int main(int argc, char *argv[])
 
    // Free the used memory.
    delete ode_solver;
+   delete fec;
 
    return 0;
 }

fem/CMakeLists.txt

@@ -68,6 +68,7 @@ set(SRCS
   fe/fe_nd.cpp
   fe/fe_nurbs.cpp
   fe/fe_pos.cpp
+  fe/fe_rbf.cpp
   fe/fe_rt.cpp
   fe/fe_ser.cpp
   fe_coll.cpp
@@ -175,6 +176,7 @@ set(HDRS
   fe/fe_nd.hpp
   fe/fe_nurbs.hpp
   fe/fe_pos.hpp
+  fe/fe_rbf.hpp
   fe/fe_rt.hpp
   fe/fe_ser.hpp
   fe_coll.hpp

fem/fe.hpp

@@ -26,5 +26,6 @@
 #include "fe/fe_nurbs.hpp"
 #include "fe/fe_pos.hpp"
 #include "fe/fe_ser.hpp"
+#include "fe/fe_rbf.hpp"
 
 #endif

fem/fe/fe_base.hpp

@@ -1026,9 +1026,6 @@ private:
 
    static Array2D<int> binom;
 
-   static void CalcMono(const int p, const double x, double *u);
-   static void CalcMono(const int p, const double x, double *u, double *d);
-
    static void CalcChebyshev(const int p, const double x, double *u);
    static void CalcChebyshev(const int p, const double x, double *u, double *d);
    static void CalcChebyshev(const int p, const double x, double *u, double *d,
@@ -1165,6 +1162,9 @@ public:
    static void CalcLegendre(const int p, const double x, double *u);
    static void CalcLegendre(const int p, const double x, double *u, double *d);
 
+   static void CalcMono(const int p, const double x, double *u);
+   static void CalcMono(const int p, const double x, double *u, double *d);
+
    ~Poly_1D();
 };
 

fem/fe/fe_rbf.hpp

@@ -0,0 +1,596 @@
+// Copyright (c) 2010-2023, 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.
+
+#ifndef MFEM_FE_RBF
+#define MFEM_FE_RBF
+
+#include "fe_base.hpp"
+
+namespace mfem
+{
+
+/** @brief Pure virtual class for dimensionless radial basis functions (RBFs).
+    Many RBFs are shaped like a Gaussian and are used here as an alternative
+    to polynomials in RBF and RK elements. The input for these
+    is a dimensionless radius r = ||x|| / h, where ||x|| is a distance and
+    h is the smoothing parameter, which controls the gradient of the RBF.
+    For reference: https://doi.org/10.1017/S0962492900000015 */
+class RBFKernel
+{
+public:
+   static const double GlobalRadius; // functions with r>=GR are considered global
+   RBFKernel() { };
+   virtual ~RBFKernel() { }
+
+   /// Evaluate the dimensionless RBF
+   virtual double BaseFunction(double r) const = 0;
+
+   /// Evaluate the derivative of the dimensionless RBF with respect to r
+   virtual double BaseDerivative(double r) const = 0;
+
+   /// Evaluate the second derivative of the dimensionless RBF with respect to r
+   virtual double BaseDerivative2(double r) const = 0;
+
+   /// The support radius, outside of which the function is zero if the function has compact support
+   virtual double Radius() const { return GlobalRadius; }
+
+   /// Does function have compact support?
+   virtual bool CompactSupport() const { return false; }
+
+   /** #brief This normalizes the smoothing parameter h such that h doesn't
+       need to be changed based on the choice of basis function */
+   virtual double HNorm() const = 0;
+};
+
+/// Gaussian RBF, exp(-r^2)
+class GaussianRBF : public RBFKernel
+{
+   // hNorm minimizes integral of Gaussian minus Wendland kernel over r=0,1
+   static const double hNorm;
+public:
+   GaussianRBF() { };
+   virtual ~GaussianRBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double HNorm() const { return hNorm; }
+};
+
+/// Multiquadric RBF, sqrt(1+r^2)
+class MultiquadricRBF : public RBFKernel
+{
+   // Same as inverse multiquadric
+   static const double hNorm;
+public:
+   MultiquadricRBF() { };
+   virtual ~MultiquadricRBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double HNorm() const { return hNorm; }
+};
+
+/// Inverse multiquadric RBF, 1/sqrt(1+r^2)
+class InvMultiquadricRBF : public RBFKernel
+{
+   // hNorm minimizes integral of Gaussian minus InvMQ kernel over r=0,0.5
+   static const double hNorm;
+public:
+   InvMultiquadricRBF() { };
+   virtual ~InvMultiquadricRBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double HNorm() const { return hNorm; }
+};
+
+/** @brief Identitcal to the Gaussian RBF, but subtracted by a factor
+    such that the function is exactly zero at the chosen radius */
+class CompactGaussianRBF : public RBFKernel
+{
+   static const double hNorm;
+   const double radius;
+   double multK, shiftK;
+
+public:
+   CompactGaussianRBF(const double rad = 5.0);
+   virtual ~CompactGaussianRBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+
+   virtual double HNorm() const { return hNorm; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+/// Identical to the Gaussian, but truncated (set to zero) at the chosen radius
+class TruncatedGaussianRBF : public RBFKernel
+{
+   static const double hNorm;
+   const double radius;
+public:
+   TruncatedGaussianRBF(const double rad = 5.0)
+      : radius(rad) { }
+   virtual ~TruncatedGaussianRBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+
+   virtual double HNorm() const { return hNorm; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+/// Wendland 11 RBF, (1-r)^3 * (1+3r) if r < 1
+class Wendland11RBF : public RBFKernel
+{
+   static const double radius;
+
+public:
+   Wendland11RBF() { }
+   virtual ~Wendland11RBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+
+   virtual double HNorm() const { return 1.0 / radius; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+/// Wendland 31 RBF, (1-r)^4 * (1+4r) if r < 1
+class Wendland31RBF : public RBFKernel
+{
+   static const double radius;
+
+public:
+   Wendland31RBF() { };
+   virtual ~Wendland31RBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+
+   virtual double HNorm() const { return 1.0 / radius; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+/// Wendland 33 RBF, (1-r)^8 * (1+8r+25r^2+32r^3) if r < 1
+class Wendland33RBF : public RBFKernel
+{
+   static const double radius;
+
+public:
+   Wendland33RBF() { };
+   virtual ~Wendland33RBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+
+   virtual double HNorm() const { return 1.0 / radius; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+/// Class for storing and creating the various RBFs
+class RBFType
+{
+public:
+   /// Represent each type of function for input/output
+   enum
+   {
+      Gaussian = 0,
+      Multiquadric = 1,
+      InvMultiquadric = 2,
+      TruncatedGaussian = 3,
+      CompactGaussian = 4,
+      Wendland11 = 5,
+      Wendland31 = 6,
+      Wendland33 = 7,
+      NumRBFTypes = 8
+   };
+
+   /// Return the requested RBF
+   static RBFKernel *GetRBF(const int rbfType)
+   {
+      switch (rbfType)
+      {
+         case RBFType::Gaussian:
+            return new GaussianRBF();
+         case RBFType::Multiquadric:
+            return new MultiquadricRBF();
+         case RBFType::InvMultiquadric:
+            return new InvMultiquadricRBF();
+         case RBFType::TruncatedGaussian:
+            return new TruncatedGaussianRBF();
+         case RBFType::CompactGaussian:
+            return new CompactGaussianRBF();
+         case RBFType::Wendland11:
+            return new Wendland11RBF();
+         case RBFType::Wendland31:
+            return new Wendland31RBF();
+         case RBFType::Wendland33:
+            return new Wendland33RBF();
+      }
+      MFEM_ABORT("unknown RBF type");
+      return NULL;
+   }
+
+   /// Abort if rbfType is invalid
+   static int Check(const int rbfType)
+   {
+      MFEM_VERIFY(0 <= rbfType && rbfType < NumRBFTypes,
+                  "unknown RBF type: " << rbfType);
+      return rbfType;
+   }
+
+   /// Convert rbf int to identifier for storage
+   static char GetChar(const int rbfType)
+   {
+      static const char ident[] = { 'G', 'M', 'I',
+                                    'T', 'C',
+                                    '1', '3', '6'
+                                  };
+      return ident[Check(rbfType)];
+   }
+
+   /// Convert identifier to rbf int
+   static int GetType(const char rbfIdent)
+   {
+      switch (rbfIdent)
+      {
+         case 'G': return Gaussian;
+         case 'M': return Multiquadric;
+         case 'I': return InvMultiquadric;
+         case 'T': return TruncatedGaussian;
+         case 'C': return CompactGaussian;
+         case '1': return Wendland11;
+         case '3': return Wendland31;
+         case '6': return Wendland33;
+      }
+      MFEM_ABORT("unknown RBF identifier: " << rbfIdent);
+      return -1;
+   }
+};
+
+/// Dimensionless distance metrics, whose output is the input for RBFs
+class DistanceMetric
+{
+protected:
+   int dim;
+public:
+   /** @brief Create a distance metric
+       @param D    Reference space dimension */
+   DistanceMetric(int D) { dim = D; }
+   virtual ~DistanceMetric() { }
+
+   /// Set the reference dimension
+   virtual void SetDim(int D) { dim = D; }
+
+   /// Given a vector of length D, return a physical distance
+   virtual void Distance(const Vector &x,
+                         double &r) const = 0;
+
+   /** @brief Given a vector of length D, return the gradient of the
+       distance with respect to the original coordinates */
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const = 0;
+
+   /** @brief Given a vector of length D, return the Hessian of the
+       distance with respect to the original coordinates */
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const = 0;
+
+   /// Create an Lp distance metric for the requested dimension and norm
+   static DistanceMetric *GetDistance(int dim, int pnorm);
+};
+
+/// Dimensionless distance with r = |x| + |y| + ...
+class L1Distance : public DistanceMetric
+{
+public:
+   L1Distance(int D) : DistanceMetric(D) { };
+   virtual ~L1Distance() { }
+
+   virtual void Distance(const Vector &x,
+                         double &r) const;
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const;
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const;
+};
+
+/// Dimensionless distance with r = (x^2 + y^2 + ...)^(1/2)
+class L2Distance : public DistanceMetric
+{
+public:
+   L2Distance(int D) : DistanceMetric(D) { };
+   virtual ~L2Distance() { }
+
+   virtual void Distance(const Vector &x,
+                         double &r) const;
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const;
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const;
+};
+
+/// Dimensionless distance with r = (x^p + y^p + ...)^(1/p)
+class LpDistance : public DistanceMetric
+{
+   const int p;
+   const double pinv;
+public:
+   LpDistance(int D, int pnorm)
+      : DistanceMetric(D),
+        p(pnorm),
+        pinv(1. / static_cast<double>(p))
+   { };
+   virtual ~LpDistance() { }
+
+   virtual void Distance(const Vector &x,
+                         double &r) const;
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const;
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const;
+};
+
+
+/** @brief Pure virtual class for a finite element with radial basis functions
+    instead of polynomials inside each element */
+class KernelFiniteElement : public ScalarFiniteElement
+{
+private:
+   // Choose whether to interpolate or project when Project is called
+   bool interpolate = false;
+public:
+   KernelFiniteElement(int D, Geometry::Type G, int Do, int O, int F)
+      : ScalarFiniteElement(D, G, Do, O, F) { }
+   virtual ~KernelFiniteElement() { }
+
+   /// Converts integration rule to vector
+   virtual void IntRuleToVec(const IntegrationPoint &ip,
+                             Vector &vec) const;
+
+   /// Is base RBF compact?
+   virtual bool IsCompact() const = 0;
+
+   /// Return base kernel
+   virtual const RBFKernel *Kernel() const = 0;
+
+   /** @brief Return whether shape function is a tensor product, used in providing indices for compact support */
+   virtual bool TensorIndexed() const { return false; }
+
+   /** @brief Get range of indices (start and end for each dimension) that
+       that are nonzero for compact support for the given integration point */
+   virtual void GetTensorIndices(const Vector &ip,
+                                 int (&indices)[3][2]) const
+   { MFEM_ABORT("GetTensorIndices(...)"); }
+
+   /** @brief Return total number of points in each dimension for tensor-
+       indexed points */
+   virtual void GetTensorNumPoints(int (&tNumPoints)[3]) const
+   { MFEM_ABORT("GetTensorNumPoints(...)"); }
+
+   using FiniteElement::Project;
+
+   virtual void Project(Coefficient &coeff, ElementTransformation &Trans,
+                        Vector &dofs) const;
+
+   virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans,
+                        Vector &dofs) const;
+
+   virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
+                        DenseMatrix &I) const;
+
+   virtual void GetLocalInterpolation(ElementTransformation &Trans,
+                                      DenseMatrix &I) const
+   { ScalarLocalInterpolation(Trans, I, *this); }
+
+   virtual void GetTransferMatrix(const FiniteElement &fe,
+                                  ElementTransformation &Trans,
+                                  DenseMatrix &I) const
+   { CheckScalarFE(fe).ScalarLocalInterpolation(Trans, I, *this); }
+};
+
+/** @brief Finite element using base radial basis functions without
+    polynomial corrections. */
+class RBFFiniteElement : public KernelFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable double r_scr, df_scr, ddf_scr;
+   mutable Vector x_scr, y_scr, dy_scr, dr_scr;
+   mutable DenseMatrix ddr_scr;
+   mutable int cInd[3][2];
+#endif
+   bool isCompact;
+   int dimPoints[3];
+   int numPointsD;
+   double delta; // Distance between points
+   double h;
+   double hPhys; // Shape parameter times distance between points times HNorm
+   double hPhysInv; // Inverse hPhys
+   double radPhys; // Radius adjusted by h
+   double faceFactor;
+   const RBFKernel *rbf;
+   const DistanceMetric *distance;
+   void InitializeGeometry();
+
+   // Get the dimensionless distance from x to the center of the RBF indexed i
+   virtual void DistanceVec(const int i,
+                            const Vector &x,
+                            Vector &y) const;
+
+public:
+   /** @brief Construct RBFFiniteElement
+       @param D           Reference space dimension
+       @param numPointsD  Number of points across the element in each dimension
+       @param rbfType     Type of radial basis function, from RBFType
+       @param distNorm    Norm used for distance, usually 2 = Euclidean distance
+       @param intOrder    Number of integration points per RBF point in each dimension
+       @param h           Shape parameter, approximately equal to the number of points in the support radius in one dimension
+       @param faceFactor  1.0 = points end on face, 0.0 = points end at dx/2 from face
+   */
+   RBFFiniteElement(const int D,
+                    const int numPointsD,
+                    const int rbfType,
+                    const int distNorm,
+                    const int intOrder,
+                    const double h,
+                    const double faceFactor);
+   virtual ~RBFFiniteElement() { delete rbf; delete distance; }
+
+   virtual bool TensorIndexed() const { return true; }
+   virtual void GetCompactIndices(const Vector &ip,
+                                  int (&indices)[3][2]) const;
+   virtual void GetGlobalIndices(const Vector &ip,
+                                 int (&indices)[3][2]) const;
+   virtual void GetTensorIndices(const Vector &ip,
+                                 int (&indices)[3][2]) const;
+   virtual void GetTensorNumPoints(int (&tNumPoints)[3]) const
+   {
+      tNumPoints[0] = dimPoints[0];
+      tNumPoints[1] = dimPoints[1];
+      tNumPoints[2] = dimPoints[2];
+   }
+
+   virtual bool IsCompact() const { return isCompact; }
+   virtual const RBFKernel *Kernel() const { return rbf; }
+
+   virtual void CalcShape(const IntegrationPoint &ip,
+                          Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &hess) const;
+};
+
+/** @brief Reproducing kernel finite element, which includes polynomial
+    corrections to the standard radial basis function finite element
+    to guarantee a chosen order of accuracy */
+class RKFiniteElement : public KernelFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable double f_scr;
+   mutable Vector x_scr, y_scr, g_scr, c_scr, s_scr, p_scr, df_scr;
+   mutable DenseMatrix q_scr, dq_scr, M_scr;
+   mutable Vector dc_scr[3], dp_scr[3];
+   mutable DenseMatrix dM_scr[3];
+   mutable DenseMatrixInverse Minv_scr;
+   mutable int cInd[3][2];
+   mutable int dimPoints[3];
+#endif
+   int polyOrd, numPoly, numPoly1d;
+   KernelFiniteElement *baseFE;
+
+   // Get the vector of polynomials for the corrections, evaluated at x
+   virtual void GetPoly(const Vector &x,
+                        Vector &p) const;
+   virtual void GetDPoly(const Vector &x,
+                         Vector &p,
+                         Vector (&dp)[3]) const;
+
+   /* Helper functions that return pieces of the RK evaluation
+      W_{RK,i} = P_i^T C_i W_{RBF,i}, where
+      P_i = [1, x, y, z, ...] is the polynomial vector evaluated at the point i,
+      M_i = P_i P_i^T W_{RBF,i} is a matrix used in calculating the corrections,
+      C_i = M_i^{-1} G are the RK corrections,
+      G = [1, 0, 0, ...] is a convencience vector */
+   virtual void GetG(Vector &g) const;
+   virtual void GetM(const Vector &baseShape,
+                     const IntegrationPoint &ip,
+                     DenseMatrix &M) const;
+   virtual void GetDM(const Vector &baseShape,
+                      const DenseMatrix &baseDeriv,
+                      const IntegrationPoint &ip,
+                      DenseMatrix &M,
+                      DenseMatrix (&dM)[3]) const;
+   virtual void AddToM(const Vector &p,
+                       const double &f,
+                       DenseMatrix &M) const;
+   virtual void AddToDM(const Vector &p,
+                        const Vector (&dp)[3],
+                        const double &f,
+                        const Vector &df,
+                        DenseMatrix (&dM)[3]) const;
+
+   // Given the corrections and base values, calculate the RK value at the ip
+   virtual void CalculateValues(const Vector &c,
+                                const Vector &baseShape,
+                                const IntegrationPoint &ip,
+                                Vector &shape) const;
+   virtual void CalculateDValues(const Vector &c,
+                                 const Vector (&dc)[3],
+                                 const Vector &baseShape,
+                                 const DenseMatrix &baseDShape,
+                                 const IntegrationPoint &ip,
+                                 DenseMatrix &dshape) const;
+
+   // The corrections make the distance dimensionless, so no smoothing parameter
+   virtual void DistanceVec(const int i,
+                            const Vector &x,
+                            Vector &y) const;
+
+public:
+   /** @brief Construct RBFFiniteElement
+       @param D           Reference space dimension
+       @param numPointsD  Number of points across the element in each dimension
+       @param rbfType     Type of radial basis function, from RBFType
+       @param distNorm    Norm used for distance, usually 2 = Euclidean distance
+       @param order       Order of polynomial correction, >= 0
+       @param intOrder    Number of integration points per RBF point in each dimension
+       @param h           Shape parameter, approximately equal to the number of points in the support radius in one dimension
+       @param faceFactor  1.0 = points end on face, 0.0 = points end at dx/2 from face
+   */
+   RKFiniteElement(const int D,
+                   const int numPointsD,
+                   const int rbfType,
+                   const int distNorm,
+                   const int order,
+                   const int intOrder,
+                   const double h,
+                   const double faceFactor);
+   virtual ~RKFiniteElement() { delete baseFE; }
+
+   virtual bool IsCompact() const { return baseFE->IsCompact(); }
+   virtual const RBFKernel *Kernel() const { return baseFE->Kernel(); }
+
+   static int GetNumPoly(int polyOrd, int dim);
+
+   virtual void CalcShape(const IntegrationPoint &ip,
+                          Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+};
+
+} // end namespace mfem
+
+#endif

fem/fe_coll.cpp

@@ -347,6 +347,36 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
          fec = new NURBSFECollection();
       }
    }
+   else if (!strncmp(name, "RBF", 3) || !strncmp(name, "RK", 2))
+   {
+      // Example: RK4_G_2_V_2D_0020_1.00_4.01
+      // (RK order 4, Gaussian, L2 dist, Value map, 2 dimensions,
+      // 20 points across element, face factor 1.0, smoothing length of 4.01)
+      const int dim = atoi(name + 10);
+      const int numPoints = atoi(name + 13);
+      const double h = atof(name + 23);
+      const double faceFactor = atof(name + 18);
+      const int rbfType = RBFType::GetType(name[4]);
+      const int distNorm = atoi(name + 6);
+      const int mapType = (name[8] == 'V'
+                           ? FiniteElement::VALUE
+                           : FiniteElement::INTEGRAL);
+      const int intOrder = 2; // keep this fixed for now
+
+      if (!strncmp(name, "RK", 2))
+      {
+         int order = atoi(name + 2);
+         fec = new LocalKernelFECollection(dim, numPoints, rbfType,
+                                           order, h, faceFactor,
+                                           intOrder, distNorm, mapType);
+      }
+      else
+      {
+         fec = new LocalKernelFECollection(dim, numPoints, rbfType,
+                                           -1, h, faceFactor,
+                                           intOrder, distNorm, mapType);
+      }
+   }
    else
    {
       MFEM_ABORT("unknown FiniteElementCollection: " << name);
@@ -3510,4 +3540,184 @@ FiniteElementCollection *NURBSFECollection::GetTraceCollection() const
    return NULL;
 }
 
+LocalKernelFECollection::LocalKernelFECollection(const int dim,
+                                                 const int numPointsD,
+                                                 const int rbfType,
+                                                 const int order,
+                                                 const double h,
+                                                 const double faceFactor,
+                                                 const int intOrder,
+                                                 const int distNorm,
+                                                 const int mapType)
+{
+   const char *mapStr = NULL;
+   switch (mapType)
+   {
+      case FiniteElement::VALUE:    mapStr = "V"; break;
+      case FiniteElement::INTEGRAL: mapStr = "I"; break;
+      default:
+         MFEM_ABORT("invalid mapType: " << mapType);
+   }
+   if (order == -1)
+   {
+      snprintf(d_name, 32, "RBF_%c_%d_%s_%dD_%04d_%1.2f_%.2f",
+               (int)RBFType::GetChar(rbfType), distNorm,
+               mapStr, dim, numPointsD, faceFactor, h);
+   }
+   else if (order >= 0)
+   {
+      snprintf(d_name, 32, "RK%d_%c_%d_%s_%dD_%04d_%1.2f_%.2f", order,
+               (int)RBFType::GetChar(rbfType), distNorm,
+               mapStr, dim, numPointsD, faceFactor, h);
+   }
+   else
+   {
+      MFEM_ABORT("invalid order: " << order);
+   }
+
+   for (int g = 0; g < Geometry::NumGeom; ++g)
+   {
+      L2_Elements[g] = NULL;
+      Tr_Elements[g] = NULL;
+   }
+   for (int i = 0; i < 2; i++)
+   {
+      SegDofOrd[i] = NULL;
+   }
+   OtherDofOrd = NULL;
+
+   if (dim == 0)
+   {
+      L2_Elements[Geometry::POINT] = new PointFiniteElement;
+   }
+   else if (dim == 1)
+   {
+      if (order == -1)
+      {
+         L2_Elements[Geometry::SEGMENT]
+            = new RBFFiniteElement(1, numPointsD, rbfType, distNorm,
+                                   intOrder, h, faceFactor);
+      }
+      else
+      {
+         L2_Elements[Geometry::SEGMENT]
+            = new RKFiniteElement(1, numPointsD, rbfType, distNorm,
+                                  order, intOrder, h, faceFactor);
+      }
+      L2_Elements[Geometry::SEGMENT]->SetMapType(mapType);
+      Tr_Elements[Geometry::POINT] = new PointFiniteElement;
+   }
+   else if (dim == 2)
+   {
+      if (order == -1)
+      {
+         L2_Elements[Geometry::SQUARE]
+            = new RBFFiniteElement(2, numPointsD, rbfType, distNorm,
+                                   intOrder, h, faceFactor);
+         Tr_Elements[Geometry::SEGMENT]
+            = new RBFFiniteElement(1, numPointsD, rbfType, distNorm,
+                                   intOrder, h, faceFactor);
+      }
+      else
+      {
+         L2_Elements[Geometry::SQUARE]
+            = new RKFiniteElement(2, numPointsD, rbfType, distNorm,
+                                  order, intOrder, h, faceFactor);
+         Tr_Elements[Geometry::SEGMENT]
+            = new RKFiniteElement(1, numPointsD, rbfType, distNorm,
+                                  order, intOrder, h, faceFactor);
+      }
+      L2_Elements[Geometry::SQUARE]->SetMapType(mapType);
+   }
+   else if (dim == 3)
+   {
+      if (order == -1)
+      {
+         L2_Elements[Geometry::CUBE]
+            = new RBFFiniteElement(3, numPointsD, rbfType, distNorm,
+                                   intOrder, h, faceFactor);
+         Tr_Elements[Geometry::SQUARE]
+            = new RBFFiniteElement(2, numPointsD, rbfType, distNorm,
+                                   intOrder, h, faceFactor);
+      }
+      else
+      {
+         L2_Elements[Geometry::CUBE]
+            = new RKFiniteElement(3, numPointsD, rbfType, distNorm,
+                                  order, intOrder, h, faceFactor);
+         Tr_Elements[Geometry::SQUARE]
+            = new RKFiniteElement(2, numPointsD, rbfType, distNorm,
+                                  order, intOrder, h, faceFactor);
+      }
+      L2_Elements[Geometry::CUBE]->SetMapType(mapType);
+   }
+
+   if (dim == 1)
+   {
+      SegDofOrd[0] = new int[2*numPointsD];
+      SegDofOrd[1] = SegDofOrd[0] + numPointsD;
+      for (int i = 0; i < numPointsD; ++i)
+      {
+         SegDofOrd[0][i] = i;
+         SegDofOrd[1][i] = numPointsD - i - 1;
+      }
+   }
+   else
+   {
+      const int geomType = TensorBasisElement::GetTensorProductGeometry(dim);
+      const int dof = L2_Elements[geomType]->GetDof();
+      OtherDofOrd = new int[dof];
+      for (int i = 0; i < dof; ++i)
+      {
+         OtherDofOrd[i] = i;
+      }
+   }
 }
+
+LocalKernelFECollection::~LocalKernelFECollection()
+{
+   delete [] OtherDofOrd;
+   delete [] SegDofOrd[0];
+   for (int i = 0; i < Geometry::NumGeom; ++i)
+   {
+      delete L2_Elements[i];
+   }
+}
+
+const FiniteElement *
+LocalKernelFECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   return L2_Elements[GeomType];
+}
+
+const FiniteElement *
+LocalKernelFECollection::TraceFiniteElementForGeometry(Geometry::Type GeomType)
+const
+{
+   return Tr_Elements[GeomType];
+}
+
+int LocalKernelFECollection::DofForGeometry(Geometry::Type GeomType) const
+{
+   if (L2_Elements[GeomType])
+   {
+      return L2_Elements[GeomType]->GetDof();
+   }
+   return 0;
+}
+
+const int *LocalKernelFECollection::DofOrderForOrientation(
+   Geometry::Type GeomType,
+   int Or) const
+{
+   if (GeomType == Geometry::SEGMENT)
+   {
+      return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
+   }
+   else
+   {
+      return (Or == 0) ? OtherDofOrd : NULL;
+   }
+}
+
+} // namespace mfem

fem/fe_coll.hpp

@@ -1305,6 +1305,64 @@ public:
    virtual int GetContType() const { return DISCONTINUOUS; }
 };
 
-}
+/** @brief Radial basis function collection.
+    Instead of using polynomials as a basis, these methods use
+    Gaussian-like functions centered at chosen points. Reproducing
+    kernels augment the radial basis functions so that the resultant
+    basis can exactly represent functions of a chosen polynomial degree. */
+class LocalKernelFECollection : public FiniteElementCollection
+{
+private:
+   int maxDim;
+   char d_name[32];
+   ScalarFiniteElement *Tr_Elements[Geometry::NumGeom];
+   ScalarFiniteElement *L2_Elements[Geometry::NumGeom];
+   int *SegDofOrd[2]; // for rotating segment dofs in 1D
+   int *OtherDofOrd;
+
+public:
+   /** @brief Create a radial basis function (RBF) or reproducing kernel (RK)
+       finite element collection.
+       @param D           Reference space dimension
+       @param numPointsD  Number of points across the element in each dimension
+       @param rbfType     Type of radial basis function, from RBFType
+       @param order       Order of polynomial correction, >= 0 for RK or -1 for
+                          RBF
+       @param h           Shape parameter, approximately equal to the number of
+                          points in the support radius in one dimension, should
+                          generally be > order
+       @param distNorm    Norm used for distance, usually 2 = Euclidean distance
+       @param intOrder    Number of integration points per RBF point in each
+                          dimension
+       @param faceFactor  1.0 = points end on face, 0.0 = points end at dx/2
+                          from face
+   */
+   LocalKernelFECollection(const int D,
+                           const int numPointsD,
+                           const int rbfType,
+                           const int order,
+                           const double h,
+                           const double faceFactor = 0.0,
+                           const int intOrder = 2,
+                           const int distNorm = 2,
+                           const int mapType = FiniteElement::VALUE);
+   virtual ~LocalKernelFECollection();
+
+   virtual const FiniteElement *
+   FiniteElementForGeometry(Geometry::Type GeomType) const;
+
+   virtual const FiniteElement *
+   TraceFiniteElementForGeometry(Geometry::Type GeomType) const;
+
+   virtual int DofForGeometry(Geometry::Type GeomType) const;
+
+   virtual const int * DofOrderForOrientation(Geometry::Type GeomType,
+                                              int Or) const;
+
+   virtual const char * Name() const { return d_name; }
+   virtual int GetContType() const { return DISCONTINUOUS; }
+};
+
+} // namespace mfem
 
 #endif