Merge branch 'complex-coef-dev' of github.com:mfem/mfem into complex-coef-dev

NURBS surface interpolation minapp

.gitignore

@@ -300,6 +300,7 @@ miniapps/nurbs/nurbs_solenoidal
 miniapps/nurbs/nurbs_printfunc
 miniapps/nurbs/nurbs_patch_ex1
 miniapps/nurbs/nurbs_curveint
+miniapps/nurbs/nurbs_surface
 miniapps/nurbs/refined.mesh
 miniapps/nurbs/mesh.*
 miniapps/nurbs/sol_?.gf
@@ -318,6 +319,7 @@ miniapps/nurbs/nurbs_naca_cmesh
 miniapps/nurbs/naca-cmesh.mesh
 miniapps/nurbs/glvis_naca-cmesh.mesh
 miniapps/nurbs/Naca_cmesh
+miniapps/nurbs/*-Surface.mesh
 
 miniapps/performance/ex1
 miniapps/performance/ex1p

CHANGELOG

@@ -29,9 +29,14 @@ Discretization improvements
 
 Meshing improvements
 --------------------
+
 - Added support for higher order meshes in Mesh::MakeSimplicial and
   ParMesh::MakeSimplicial.
 
+- Added a new miniapp for interpolating a surface grid of points in 3D using a
+  smooth NURBS surface, that can then be sampled at arbitrary resolution while
+  staying close to the original geometry. See miniapps/nurbs/nurbs_surface.
+
 GPU computing
 -------------
 - The function Vector::SetSubVector(const Array<int> &, const real_t) now

examples/ex22.cpp

@@ -62,14 +62,9 @@ static real_t epsilon_ = 1.0;
 static real_t sigma_ = 20.0;
 static real_t omega_ = 10.0;
 
-real_t u0_real_exact(const Vector &);
-real_t u0_imag_exact(const Vector &);
-
-void u1_real_exact(const Vector &, Vector &);
-void u1_imag_exact(const Vector &, Vector &);
-
-void u2_real_exact(const Vector &, Vector &);
-void u2_imag_exact(const Vector &, Vector &);
+complex<real_t> u0_exact(const Vector &x);
+void u1_exact(const Vector &, ComplexVector &);
+void u2_exact(const Vector &, ComplexVector &);
 
 bool check_for_inline_mesh(const char * mesh_file);
 
@@ -215,54 +210,48 @@ int main(int argc, char *argv[])
    ComplexGridFunction * u_exact = NULL;
    if (exact_sol) { u_exact = new ComplexGridFunction(fespace); }
 
-   FunctionCoefficient u0_r(u0_real_exact);
-   FunctionCoefficient u0_i(u0_imag_exact);
-   VectorFunctionCoefficient u1_r(dim, u1_real_exact);
-   VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
-   VectorFunctionCoefficient u2_r(dim, u2_real_exact);
-   VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
+   ComplexFunctionCoefficient u0(u0_exact);
+   ComplexVectorFunctionCoefficient u1(dim, u1_exact);
+   ComplexVectorFunctionCoefficient u2(dim, u2_exact);
 
-   ConstantCoefficient zeroCoef(0.0);
-   ConstantCoefficient oneCoef(1.0);
+   ComplexConstantCoefficient oneCoef(1.0);
 
-   Vector zeroVec(dim); zeroVec = 0.0;
    Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
-   VectorConstantCoefficient zeroVecCoef(zeroVec);
-   VectorConstantCoefficient oneVecCoef(oneVec);
+   ComplexVectorConstantCoefficient oneVecCoef(oneVec);
 
    switch (prob)
    {
       case 0:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficient(u0_r, u0_i, ess_bdr);
-            u_exact->ProjectCoefficient(u0_r, u0_i);
+            u.ProjectBdrCoefficient(u0, ess_bdr);
+            u_exact->ProjectCoefficient(u0);
          }
          else
          {
-            u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
+            u.ProjectBdrCoefficient(oneCoef, ess_bdr);
          }
          break;
       case 1:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientTangent(u1_r, u1_i, ess_bdr);
-            u_exact->ProjectCoefficient(u1_r, u1_i);
+            u.ProjectBdrCoefficientTangent(u1, ess_bdr);
+            u_exact->ProjectCoefficient(u1);
          }
          else
          {
-            u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientTangent(oneVecCoef, ess_bdr);
          }
          break;
       case 2:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientNormal(u2_r, u2_i, ess_bdr);
-            u_exact->ProjectCoefficient(u2_r, u2_i);
+            u.ProjectBdrCoefficientNormal(u2, ess_bdr);
+            u_exact->ProjectCoefficient(u2);
          }
          else
          {
-            u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientNormal(oneVecCoef, ess_bdr);
          }
          break;
       default: break; // This should be unreachable
@@ -300,27 +289,24 @@ int main(int argc, char *argv[])
    ConstantCoefficient lossCoef(omega_ * sigma_);
    ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
 
+   ComplexConstantCoefficient complexMassCoef(-omega_ * omega_ * epsilon_,
+                                              omega_ * sigma_);
+
    SesquilinearForm *a = new SesquilinearForm(fespace, conv);
    if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    switch (prob)
    {
       case 0:
-         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new MassIntegrator(massCoef),
-                                new MassIntegrator(lossCoef));
+         a->AddDomainIntegrator<DiffusionIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<MassIntegrator>(complexMassCoef);
          break;
       case 1:
-         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
+         a->AddDomainIntegrator<CurlCurlIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
          break;
       case 2:
-         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
+         a->AddDomainIntegrator<DivDivIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
          break;
       default: break; // This should be unreachable
    }
@@ -436,29 +422,24 @@ int main(int argc, char *argv[])
 
    if (exact_sol)
    {
-      real_t err_r = -1.0;
-      real_t err_i = -1.0;
+      real_t err_u = -1.0;
 
       switch (prob)
       {
          case 0:
-            err_r = u.real().ComputeL2Error(u0_r);
-            err_i = u.imag().ComputeL2Error(u0_i);
+            err_u = u.ComputeL2Error(u0);
             break;
          case 1:
-            err_r = u.real().ComputeL2Error(u1_r);
-            err_i = u.imag().ComputeL2Error(u1_i);
+            err_u = u.ComputeL2Error(u1);
             break;
          case 2:
-            err_r = u.real().ComputeL2Error(u2_r);
-            err_i = u.imag().ComputeL2Error(u2_i);
+            err_u = u.ComputeL2Error(u2);
             break;
          default: break; // This should be unreachable
       }
 
       cout << endl;
-      cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
-      cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
+      cout << "|| u_h - u ||_{L^2} = " << err_u << endl;
       cout << endl;
    }
 
@@ -564,36 +545,14 @@ complex<real_t> u0_exact(const Vector &x)
    return std::exp(-i * kappa * x[dim - 1]);
 }
 
-real_t u0_real_exact(const Vector &x)
-{
-   return u0_exact(x).real();
-}
-
-real_t u0_imag_exact(const Vector &x)
-{
-   return u0_exact(x).imag();
-}
-
-void u1_real_exact(const Vector &x, Vector &v)
+void u1_exact(const Vector &x, ComplexVector &v)
 {
    int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
+   v.SetSize(dim); v = 0.0; v[0] = u0_exact(x);
 }
 
-void u1_imag_exact(const Vector &x, Vector &v)
+void u2_exact(const Vector &x, ComplexVector &v)
 {
    int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
-}
-
-void u2_real_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
-}
-
-void u2_imag_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_exact(x);
 }

examples/ex22p.cpp

@@ -62,6 +62,10 @@ static real_t epsilon_ = 1.0;
 static real_t sigma_ = 20.0;
 static real_t omega_ = 10.0;
 
+complex<real_t> u0_exact(const Vector &x);
+void u1_exact(const Vector &, ComplexVector &);
+void u2_exact(const Vector &, ComplexVector &);
+
 real_t u0_real_exact(const Vector &);
 real_t u0_imag_exact(const Vector &);
 
@@ -244,13 +248,22 @@ int main(int argc, char *argv[])
    ParComplexGridFunction * u_exact = NULL;
    if (exact_sol) { u_exact = new ParComplexGridFunction(fespace); }
 
+   ComplexFunctionCoefficient u0(u0_exact);
+   ComplexVectorFunctionCoefficient u1(dim, u1_exact);
+   ComplexVectorFunctionCoefficient u2(dim, u2_exact);
+
+   ComplexConstantCoefficient oneCoef(1.0);
+
+   Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
+   ComplexVectorConstantCoefficient oneVecCoef(oneVec);
+
    FunctionCoefficient u0_r(u0_real_exact);
    FunctionCoefficient u0_i(u0_imag_exact);
    VectorFunctionCoefficient u1_r(dim, u1_real_exact);
    VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
    VectorFunctionCoefficient u2_r(dim, u2_real_exact);
    VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
-
+   /*
    ConstantCoefficient zeroCoef(0.0);
    ConstantCoefficient oneCoef(1.0);
 
@@ -258,40 +271,40 @@ int main(int argc, char *argv[])
    Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
    VectorConstantCoefficient zeroVecCoef(zeroVec);
    VectorConstantCoefficient oneVecCoef(oneVec);
-
+   */
    switch (prob)
    {
       case 0:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficient(u0_r, u0_i, ess_bdr);
-            u_exact->ProjectCoefficient(u0_r, u0_i);
+            u.ProjectBdrCoefficient(u0, ess_bdr);
+            u_exact->ProjectCoefficient(u0);
          }
          else
          {
-            u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
+            u.ProjectBdrCoefficient(oneCoef, ess_bdr);
          }
          break;
       case 1:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientTangent(u1_r, u1_i, ess_bdr);
-            u_exact->ProjectCoefficient(u1_r, u1_i);
+            u.ProjectBdrCoefficientTangent(u1, ess_bdr);
+            u_exact->ProjectCoefficient(u1);
          }
          else
          {
-            u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientTangent(oneVecCoef, ess_bdr);
          }
          break;
       case 2:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientNormal(u2_r, u2_i, ess_bdr);
-            u_exact->ProjectCoefficient(u2_r, u2_i);
+            u.ProjectBdrCoefficientNormal(u2, ess_bdr);
+            u_exact->ProjectCoefficient(u2);
          }
          else
          {
-            u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientNormal(oneVecCoef, ess_bdr);
          }
          break;
       default: break; // This should be unreachable
@@ -331,27 +344,24 @@ int main(int argc, char *argv[])
    ConstantCoefficient lossCoef(omega_ * sigma_);
    ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
 
+   ComplexConstantCoefficient complexMassCoef(-omega_ * omega_ * epsilon_,
+                                              omega_ * sigma_);
+
    ParSesquilinearForm *a = new ParSesquilinearForm(fespace, conv);
    if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    switch (prob)
    {
       case 0:
-         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new MassIntegrator(massCoef),
-                                new MassIntegrator(lossCoef));
+         a->AddDomainIntegrator<DiffusionIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<MassIntegrator>(complexMassCoef);
          break;
       case 1:
-         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
+         a->AddDomainIntegrator<CurlCurlIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
          break;
       case 2:
-         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
+         a->AddDomainIntegrator<DivDivIntegrator>(stiffnessCoef);
+         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
          break;
       default: break; // This should be unreachable
    }
@@ -475,22 +485,18 @@ int main(int argc, char *argv[])
 
    if (exact_sol)
    {
-      real_t err_r = -1.0;
-      real_t err_i = -1.0;
+      real_t err_u = -1.0;
 
       switch (prob)
       {
          case 0:
-            err_r = u.real().ComputeL2Error(u0_r);
-            err_i = u.imag().ComputeL2Error(u0_i);
+            err_u = u.ComputeL2Error(u0);
             break;
          case 1:
-            err_r = u.real().ComputeL2Error(u1_r);
-            err_i = u.imag().ComputeL2Error(u1_i);
+            err_u = u.ComputeL2Error(u1);
             break;
          case 2:
-            err_r = u.real().ComputeL2Error(u2_r);
-            err_i = u.imag().ComputeL2Error(u2_i);
+            err_u = u.ComputeL2Error(u2);
             break;
          default: break; // This should be unreachable
       }
@@ -498,8 +504,7 @@ int main(int argc, char *argv[])
       if ( myid == 0 )
       {
          cout << endl;
-         cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
-         cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
+         cout << "|| u_h - u ||_{L^2} = " << err_u << endl;
          cout << endl;
       }
    }
@@ -627,6 +632,12 @@ real_t u0_imag_exact(const Vector &x)
    return u0_exact(x).imag();
 }
 
+void u1_exact(const Vector &x, ComplexVector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_exact(x);
+}
+
 void u1_real_exact(const Vector &x, Vector &v)
 {
    int dim = x.Size();
@@ -639,6 +650,12 @@ void u1_imag_exact(const Vector &x, Vector &v)
    v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
 }
 
+void u2_exact(const Vector &x, ComplexVector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_exact(x);
+}
+
 void u2_real_exact(const Vector &x, Vector &v)
 {
    int dim = x.Size();

fem/CMakeLists.txt

@@ -59,6 +59,7 @@ set(SRCS
   integ/nonlininteg_vecconvection_pa.cpp
   integ/nonlininteg_vecconvection_mf.cpp
   coefficient.cpp
+  complex_coefficient.cpp
   complex_fem.cpp
   convergence.cpp
   datacollection.cpp
@@ -176,6 +177,7 @@ set(HDRS
   integ/bilininteg_hcurlhdiv_kernels.hpp
   integ/bilininteg_mass_kernels.hpp
   coefficient.hpp
+  complex_coefficient.hpp
   complex_fem.hpp
   convergence.hpp
   datacollection.hpp

fem/complex_coefficient.cpp

@@ -0,0 +1,217 @@
+// Copyright (c) 2010-2025, 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 "complex_fem.hpp"
+#include "../general/forall.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+real_t
+RealPartCoefficient::Eval(ElementTransformation &T,
+                          const IntegrationPoint &ip)
+{
+   complex_t val = complex_coef_.Eval(T, ip);
+   return val.real();
+}
+
+real_t
+ImagPartCoefficient::Eval(ElementTransformation &T,
+                          const IntegrationPoint &ip)
+{
+   complex_t val = complex_coef_.Eval(T, ip);
+   return val.imag();
+}
+
+RealPartVectorCoefficient::RealPartVectorCoefficient(ComplexVectorCoefficient &
+                                                     complex_vcoef)
+   : VectorCoefficient(complex_vcoef.GetVDim()),
+     complex_vcoef_(complex_vcoef),
+     val_(vdim)
+{}
+
+void
+RealPartVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
+                                const IntegrationPoint &ip)
+{
+   complex_vcoef_.Eval(val_, T, ip);
+   V = val_.real();
+}
+
+ImagPartVectorCoefficient::ImagPartVectorCoefficient(ComplexVectorCoefficient &
+                                                     complex_vcoef)
+   : VectorCoefficient(complex_vcoef.GetVDim()),
+     complex_vcoef_(complex_vcoef),
+     val_(vdim)
+{}
+
+void
+ImagPartVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
+                                const IntegrationPoint &ip)
+{
+   complex_vcoef_.Eval(val_, T, ip);
+   V = val_.imag();
+}
+
+RealPartMatrixCoefficient::RealPartMatrixCoefficient(ComplexMatrixCoefficient &
+                                                     complex_mcoef)
+   : MatrixCoefficient(complex_mcoef.GetHeight(), complex_mcoef.GetWidth()),
+     complex_mcoef_(complex_mcoef),
+     val_(height, width)
+{}
+
+void
+RealPartMatrixCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
+                                const IntegrationPoint &ip)
+{
+   complex_mcoef_.Eval(val_, T, ip);
+   M = val_.real();
+}
+
+ImagPartMatrixCoefficient::ImagPartMatrixCoefficient(ComplexMatrixCoefficient &
+                                                     complex_mcoef)
+   : MatrixCoefficient(complex_mcoef.GetHeight(), complex_mcoef.GetWidth()),
+     complex_mcoef_(complex_mcoef),
+     val_(height, width)
+{}
+
+void
+ImagPartMatrixCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
+                                const IntegrationPoint &ip)
+{
+   complex_mcoef_.Eval(val_, T, ip);
+   M = val_.imag();
+}
+
+ComplexCoefficient::ComplexCoefficient()
+   : time(0.),
+     re_part_coef_(*this), im_part_coef_(*this),
+     real_coef_(re_part_coef_), imag_coef_(im_part_coef_)
+{ }
+
+ComplexCoefficient::ComplexCoefficient(Coefficient &c_r,
+                                       Coefficient &c_i)
+   : time(c_r.GetTime()),
+     re_part_coef_(*this), im_part_coef_(*this),
+     real_coef_(c_r), imag_coef_(c_i)
+{
+   c_i.SetTime(time);
+}
+
+complex_t
+ComplexCoefficient::Eval(ElementTransformation &T,
+                         const IntegrationPoint &ip)
+{
+   // Avoid circular dependency
+   MFEM_VERIFY(std::addressof(real_coef_) != std::addressof(re_part_coef_) &&
+               std::addressof(imag_coef_) != std::addressof(im_part_coef_),
+               "Classes dervied from ComplexCoefficient must either "
+               "implement an Eval method or supply Coefficients "
+               "for both the real and imaginary parts of the field.");
+
+   return complex_t(real_coef_.Eval(T, ip), imag_coef_.Eval(T, ip));
+}
+
+ComplexVectorCoefficient::ComplexVectorCoefficient(VectorCoefficient &v_r,
+                                                   VectorCoefficient &v_i)
+   : vdim(v_r.GetVDim()), time(v_r.GetTime()),
+     re_part_vcoef_(*this), im_part_vcoef_(*this),
+     real_vcoef_(v_r), imag_vcoef_(v_i)
+{
+   MFEM_ASSERT(v_r.GetVDim() == v_i.GetVDim(), "ComplexVectorCoefficient"
+               " - incompatible vector dimensions of real and imaginary parts.");
+
+   v_i.SetTime(time);
+}
+
+void ComplexVectorCoefficient::Eval(ComplexVector &V, ElementTransformation &T,
+                                    const IntegrationPoint &ip)
+{
+   // Avoid circular dependency
+   MFEM_VERIFY(std::addressof(real_vcoef_) != std::addressof(re_part_vcoef_) &&
+               std::addressof(imag_vcoef_) != std::addressof(im_part_vcoef_),
+               "Classes dervied from ComplexVectorCoefficient must either "
+               "implement an Eval method or supply VectorCoefficients "
+               "for both the real and imaginary parts of the field.");
+
+   V_r_.SetSize(vdim);
+   V_i_.SetSize(vdim);
+
+   real_vcoef_.Eval(V_r_, T, ip);
+   imag_vcoef_.Eval(V_i_, T, ip);
+
+   V.Set(V_r_, V_i_);
+}
+
+ComplexConstantCoefficient::ComplexConstantCoefficient(
+   const complex_t z)
+   : val(z), real_coef(z.real()), imag_coef(z.imag())
+{
+   real_coef_ = real_coef;
+   imag_coef_ = imag_coef;
+}
+
+ComplexConstantCoefficient::ComplexConstantCoefficient(
+   real_t z_r, real_t z_i)
+   : real_coef(z_r), imag_coef(z_i)
+{
+   val = complex_t(z_r, z_i);
+
+   real_coef_ = real_coef;
+   imag_coef_ = imag_coef;
+}
+
+complex_t ComplexFunctionCoefficient::Eval(ElementTransformation & T,
+                                           const IntegrationPoint & ip)
+{
+   real_t x[3];
+   Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+
+   if (Function)
+   {
+      return Function(transip);
+   }
+   else
+   {
+      return TDFunction(transip, GetTime());
+   }
+}
+
+void ComplexVectorFunctionCoefficient::Eval(ComplexVector &V,
+                                            ElementTransformation &T,
+                                            const IntegrationPoint &ip)
+{
+   real_t x[3];
+   Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+
+   V.SetSize(vdim);
+   if (Function)
+   {
+      Function(transip, V);
+   }
+   else
+   {
+      TDFunction(transip, GetTime(), V);
+   }
+   if (Q)
+   {
+      V *= Q->Eval(T, ip, GetTime());
+   }
+}
+
+} // end namespace mfem
+

fem/complex_coefficient.hpp

@@ -0,0 +1,523 @@
+// Copyright (c) 2010-2025, 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_COMPLEX_COEFFICIENT
+#define MFEM_COMPLEX_COEFFICIENT
+
+#include "../config/config.hpp"
+#include "../linalg/linalg.hpp"
+#include "coefficient.hpp"
+#include "intrules.hpp"
+#include "eltrans.hpp"
+
+namespace mfem
+{
+
+class ComplexCoefficient;
+class ComplexVectorCoefficient;
+class ComplexMatrixCoefficient;
+
+/// Standard Coefficient which returns the real part of a ComplexCoefficient
+class RealPartCoefficient : public Coefficient
+{
+private:
+   ComplexCoefficient &complex_coef_;
+
+public:
+   RealPartCoefficient(ComplexCoefficient & complex_coef)
+      : complex_coef_(complex_coef) {}
+
+   real_t Eval(ElementTransformation &T,
+               const IntegrationPoint &ip);
+};
+
+/// Standard Coefficient which returns the imaginary part of a
+/// ComplexCoefficient
+class ImagPartCoefficient : public Coefficient
+{
+private:
+   ComplexCoefficient &complex_coef_;
+
+public:
+   ImagPartCoefficient(ComplexCoefficient & complex_coef)
+      : complex_coef_(complex_coef) {}
+
+   real_t Eval(ElementTransformation &T,
+               const IntegrationPoint &ip);
+};
+
+typedef ImagPartCoefficient ImaginaryPartCoefficient;
+
+class RealPartVectorCoefficient : public VectorCoefficient
+{
+private:
+   ComplexVectorCoefficient &complex_vcoef_;
+   mutable ComplexVector val_;
+
+public:
+   RealPartVectorCoefficient(ComplexVectorCoefficient & complex_vcoef);
+
+   void Eval(Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip);
+};
+
+class ImagPartVectorCoefficient : public VectorCoefficient
+{
+private:
+   ComplexVectorCoefficient &complex_vcoef_;
+   mutable ComplexVector val_;
+
+public:
+   ImagPartVectorCoefficient(ComplexVectorCoefficient & complex_vcoef);
+
+   void Eval(Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip);
+};
+
+typedef ImagPartVectorCoefficient ImaginaryPartVectorCoefficient;
+
+class RealPartMatrixCoefficient : public MatrixCoefficient
+{
+private:
+   ComplexMatrixCoefficient &complex_mcoef_;
+   mutable ComplexTypeDenseMatrix val_;
+
+public:
+   RealPartMatrixCoefficient(ComplexMatrixCoefficient & complex_mcoef);
+
+   void Eval(DenseMatrix &M, ElementTransformation &T,
+             const IntegrationPoint &ip);
+};
+
+class ImagPartMatrixCoefficient : public MatrixCoefficient
+{
+private:
+   ComplexMatrixCoefficient &complex_mcoef_;
+   mutable ComplexTypeDenseMatrix val_;
+
+public:
+   ImagPartMatrixCoefficient(ComplexMatrixCoefficient & complex_mcoef);
+
+   void Eval(DenseMatrix &V, ElementTransformation &T,
+             const IntegrationPoint &ip);
+};
+
+typedef ImagPartMatrixCoefficient ImaginaryPartMatrixCoefficient;
+
+/** @brief Base class ComplexCoefficients that optionally depend on space and
+    time. These are used by the SesquilinearForm, ComplexLinearForm, and
+    ComplexGridFunction classes to represent the physical coefficients in
+    the PDEs that are being discretized. This class can also be used in a more
+    general way to represent functions that don't necessarily belong to a FE
+    space, e.g., to project onto ComplexGridFunctions to use as initial
+    conditions, exact solutions, etc. See, e.g., ex22 for these uses. */
+class ComplexCoefficient
+{
+protected:
+   real_t time;
+
+private:
+   RealPartCoefficient re_part_coef_;
+   ImagPartCoefficient im_part_coef_;
+
+protected:
+   Coefficient &real_coef_;
+   Coefficient &imag_coef_;
+
+public:
+
+   ComplexCoefficient();
+   ComplexCoefficient(Coefficient &c_r, Coefficient &c_i);
+
+   /// Set the time for time dependent coefficients
+   virtual void SetTime(real_t t)
+   { time = t; real_coef_.SetTime(t); imag_coef_.SetTime(t); }
+
+   /// Get the time for time dependent coefficients
+   real_t GetTime() { return time; }
+
+   /** @brief Evaluate the coefficient in the element described by @a T at the
+       point @a ip. */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   virtual complex_t Eval(ElementTransformation &T,
+                          const IntegrationPoint &ip);
+
+   /** @brief Evaluate the coefficient in the element described by @a T at the
+       point @a ip at time @a t. */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   complex_t Eval(ElementTransformation &T,
+                  const IntegrationPoint &ip, real_t t)
+   {
+      SetTime(t);
+      return Eval(T, ip);
+   }
+
+   /** @brief Access a standard Coefficient object reproducing the real part of
+       the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its real part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual Coefficient & real() { return real_coef_; }
+
+   /** @brief Access a standard Coefficient object reproducing the imaginary
+       part of the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its imaginary part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual Coefficient & imag() { return imag_coef_; }
+
+   virtual ~ComplexCoefficient() { }
+};
+
+/** @brief Base class ComplexVectorCoefficients that optionally depend
+    on space and time. These are used by the SesquilinearForm,
+    ComplexLinearForm, and ComplexGridFunction classes to represent
+    the physical vector-valued coefficients in the PDEs that are being
+    discretized. This class can also be used in a more general way to
+    represent functions that don't necessarily belong to a FE space,
+    e.g., to project onto ComplexGridFunctions to use as initial
+    conditions, exact solutions, etc. See, e.g., ex22 for these
+    uses. */
+class ComplexVectorCoefficient
+{
+protected:
+   int vdim;
+   real_t time;
+
+private:
+   RealPartVectorCoefficient re_part_vcoef_;
+   ImagPartVectorCoefficient im_part_vcoef_;
+
+protected:
+   VectorCoefficient &real_vcoef_;
+   VectorCoefficient &imag_vcoef_;
+
+   mutable Vector V_r_;
+   mutable Vector V_i_;
+
+public:
+   ComplexVectorCoefficient(int vd)
+      : vdim(vd), time(0.),
+        re_part_vcoef_(*this), im_part_vcoef_(*this),
+        real_vcoef_(re_part_vcoef_), imag_vcoef_(im_part_vcoef_)
+   { }
+
+   ComplexVectorCoefficient(VectorCoefficient &v_r, VectorCoefficient &v_i);
+
+
+   /// Set the time for time dependent coefficients
+   virtual void SetTime(real_t t)
+   { time = t; real_vcoef_.SetTime(t); imag_vcoef_.SetTime(t); }
+
+   /// Get the time for time dependent coefficients
+   real_t GetTime() { return time; }
+
+   /// Returns dimension of the vector.
+   int GetVDim() { return vdim; }
+
+   /** @brief Evaluate the vector coefficient in the element described by @a T
+       at the point @a ip, storing the result in @a V. */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   virtual void Eval(ComplexVector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+
+   /** @brief Evaluate the vector coefficient in the element described by @a T
+       at the point @a ip at time @a t, storing the result in @a V. */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   void Eval(ComplexVector &V, ElementTransformation &T,
+             const IntegrationPoint &ip, real_t t)
+   {
+      SetTime(t);
+      Eval(V, T, ip);
+   }
+
+   /** @brief Access a standard Coefficient object reproducing the real part of
+       the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its real part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual VectorCoefficient & real() { return real_vcoef_; }
+
+   /** @brief Access a standard Coefficient object reproducing the imaginary
+       part of the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its imaginary part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual VectorCoefficient & imag() { return imag_vcoef_; }
+
+   virtual ~ComplexVectorCoefficient() { }
+};
+
+/** @brief Base class ComplexMatrixCoefficients that optionally depend
+    on space and time. These are used by the SesquilinearForm,
+    ComplexLinearForm, and ComplexGridFunction classes to represent
+    the physical matrix-valued coefficients in the PDEs that are being
+    discretized. This class can also be used in a more general way to
+    represent functions that don't necessarily belong to a FE space.
+    See, e.g., ex22 for these uses. */
+class ComplexMatrixCoefficient
+{
+protected:
+   int height, width;
+   real_t time;
+
+private:
+   RealPartMatrixCoefficient re_part_mcoef_;
+   ImagPartMatrixCoefficient im_part_mcoef_;
+
+protected:
+   MatrixCoefficient &real_mcoef_;
+   MatrixCoefficient &imag_mcoef_;
+
+   mutable DenseMatrix M_r_;
+   mutable DenseMatrix M_i_;
+
+public:
+   /// Construct a dim x dim matrix coefficient.
+   explicit ComplexMatrixCoefficient(int dim)
+      : height(dim), width(dim), time(0.),
+        re_part_mcoef_(*this), im_part_mcoef_(*this),
+        real_mcoef_(re_part_mcoef_), imag_mcoef_(im_part_mcoef_)
+   { }
+
+   /// Construct a h x w matrix coefficient.
+   ComplexMatrixCoefficient(int h, int w) :
+      height(h), width(w), time(0.),
+      re_part_mcoef_(*this), im_part_mcoef_(*this),
+      real_mcoef_(re_part_mcoef_), imag_mcoef_(im_part_mcoef_)
+   { }
+
+   /// Set the time for time dependent coefficients
+   virtual void SetTime(real_t t) { time = t; }
+
+   /// Get the time for time dependent coefficients
+   real_t GetTime() { return time; }
+
+   /// Get the height of the matrix.
+   int GetHeight() const { return height; }
+
+   /// Get the width of the matrix.
+   int GetWidth() const { return width; }
+
+   /// For backward compatibility get the width of the matrix.
+   int GetVDim() const { return width; }
+
+   /** @brief Evaluate the matrix coefficient in the element described by @a T
+       at the point @a ip, storing the result in @a K. */
+   /** @note When this method is called, the caller must make sure that the
+       IntegrationPoint associated with @a T is the same as @a ip. This can be
+       achieved by calling T.SetIntPoint(&ip). */
+   virtual void Eval(ComplexTypeDenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip) = 0;
+
+   /** @brief Access a standard Coefficient object reproducing the real part of
+       the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its real part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual MatrixCoefficient & real() { return real_mcoef_; }
+
+   /** @brief Access a standard Coefficient object reproducing the imaginary
+       part of the complex-valued field */
+   /** @note By default this method returns an internal object which
+       computes the complex value using the above Eval method and
+       returns its imaginary part. Custom implementations may choose to
+       override this method with a more efficient real-valued
+       coefficient. */
+   virtual MatrixCoefficient & imag() { return imag_mcoef_; }
+
+   virtual ~ComplexMatrixCoefficient() { }
+};
+
+/// A complex-valued coefficient that is constant across space and time
+class ComplexConstantCoefficient : public ComplexCoefficient
+{
+private:
+   complex_t val;
+
+   ConstantCoefficient real_coef;
+   ConstantCoefficient imag_coef;
+
+public:
+   ComplexConstantCoefficient(const complex_t z);
+
+   ComplexConstantCoefficient(real_t z_r, real_t z_i = 0.);
+
+   complex_t Eval(ElementTransformation &T,
+                  const IntegrationPoint &ip) { return val; }
+};
+
+/// Complex-valued vector coefficient that is constant in space and time.
+class ComplexVectorConstantCoefficient : public ComplexVectorCoefficient
+{
+private:
+   ComplexVector vec;
+
+public:
+   /// Construct the coefficient with constant vector @a v.
+   ComplexVectorConstantCoefficient(const ComplexVector &v)
+      : ComplexVectorCoefficient(v.Size()), vec(v) { }
+
+   /// Construct the coefficient with constant vector @a v.
+   ComplexVectorConstantCoefficient(const Vector &v)
+      : ComplexVectorCoefficient(v.Size()), vec(v) { }
+
+   using ComplexVectorCoefficient::Eval;
+
+   ///  Evaluate the vector coefficient at @a ip.
+   void Eval(ComplexVector &V, ElementTransformation &T,
+             const IntegrationPoint &ip) override { V = vec; }
+
+   /// Return a reference to the constant vector in this class.
+   const ComplexVector& GetVec() const { return vec; }
+};
+
+/// Complex-valued vector coefficient that is constant in space and time.
+class ComplexMatrixConstantCoefficient : public ComplexMatrixCoefficient
+{
+private:
+   ComplexTypeDenseMatrix mat;
+
+public:
+   /// Construct the coefficient with constant vector @a v.
+   ComplexMatrixConstantCoefficient(const ComplexTypeDenseMatrix &m)
+      : ComplexMatrixCoefficient(m.Height(), m.Width()), mat(m) { }
+
+   /// Construct the coefficient with constant vector @a v.
+   ComplexMatrixConstantCoefficient(const DenseMatrix &m)
+      : ComplexMatrixCoefficient(m.Height(), m.Width()), mat(m) { }
+
+   using ComplexMatrixCoefficient::Eval;
+
+   ///  Evaluate the matrix coefficient at @a ip.
+   void Eval(ComplexTypeDenseMatrix &M, ElementTransformation &T,
+             const IntegrationPoint &ip) override { M = mat; }
+
+   /// Return a reference to the constant matrix in this class.
+   const ComplexTypeDenseMatrix& GetMat() const { return mat; }
+};
+
+/// A general complex-valued function coefficient
+class ComplexFunctionCoefficient : public ComplexCoefficient
+{
+protected:
+   std::function<complex_t(const Vector &)> Function;
+   std::function<complex_t(const Vector &, real_t)> TDFunction;
+
+public:
+   /// Define a time-independent coefficient from a std function
+   /** \param F time-independent std::function */
+   ComplexFunctionCoefficient(std::function<complex_t
+                              (const Vector &)> F)
+      : Function(std::move(F))
+   { }
+
+   /// Define a time-dependent coefficient from a std function
+   /** \param TDF time-dependent function */
+   ComplexFunctionCoefficient(std::function<complex_t
+                              (const Vector &, real_t)> TDF)
+      : TDFunction(std::move(TDF))
+   { }
+
+   /// (DEPRECATED) Define a time-independent coefficient from a C-function
+   /** @deprecated Use the method where the C-function, @a f, uses a const
+       Vector argument instead of Vector. */
+   MFEM_DEPRECATED ComplexFunctionCoefficient(complex_t
+                                              (*f)(Vector &))
+   {
+      // Cast first to (void*) to suppress a warning from newer version of
+      // Clang when using -Wextra.
+      Function = reinterpret_cast<complex_t(*)
+                 (const Vector&)>((void*)f);
+      TDFunction = NULL;
+   }
+
+   /// (DEPRECATED) Define a time-dependent coefficient from a C-function
+   /** @deprecated Use the method where the C-function, @a tdf, uses a const
+       Vector argument instead of Vector. */
+   MFEM_DEPRECATED ComplexFunctionCoefficient(complex_t
+                                              (*tdf)(Vector &, real_t))
+   {
+      Function = NULL;
+      // Cast first to (void*) to suppress a warning from newer version of
+      // Clang when using -Wextra.
+      TDFunction =
+         reinterpret_cast<complex_t(*)(const Vector&,
+                                       real_t)>((void*)tdf);
+   }
+
+   /// Evaluate the coefficient at @a ip.
+   complex_t Eval(ElementTransformation &T,
+                  const IntegrationPoint &ip) override;
+};
+
+/// A general vector function coefficient
+class ComplexVectorFunctionCoefficient : public ComplexVectorCoefficient
+{
+private:
+   std::function<void(const Vector &, ComplexVector &)> Function;
+   std::function<void(const Vector &, real_t, ComplexVector &)> TDFunction;
+   ComplexCoefficient *Q;
+
+public:
+   /// Define a time-independent complex-valued vector coefficient
+   /// from a std function
+   /** \param dim - the size of the vector
+       \param F - time-independent function
+       \param q - optional scalar Coefficient to scale the vector coefficient */
+   ComplexVectorFunctionCoefficient(int dim,
+                                    std::function<void(const Vector &,
+                                                       ComplexVector &)> F,
+                                    ComplexCoefficient *q = nullptr)
+      : ComplexVectorCoefficient(dim), Function(std::move(F)), Q(q)
+   { }
+
+   /// Define a time-dependent complex-valued vector coefficient from
+   /// a std function
+   /** \param dim - the size of the vector
+       \param TDF - time-dependent function
+       \param q - optional scalar ComplexCoefficient to scale the vector coefficient */
+   ComplexVectorFunctionCoefficient(int dim,
+                                    std::function<void(const Vector &, real_t,
+                                                       ComplexVector &)> TDF,
+                                    ComplexCoefficient *q = nullptr)
+      : ComplexVectorCoefficient(dim), TDFunction(std::move(TDF)), Q(q)
+   { }
+
+   using ComplexVectorCoefficient::Eval;
+   /// Evaluate the vector coefficient at @a ip.
+   void Eval(ComplexVector &V, ElementTransformation &T,
+             const IntegrationPoint &ip) override;
+
+   virtual ~ComplexVectorFunctionCoefficient() { }
+};
+
+} // end namespace mfem
+
+#endif

fem/complex_fem.cpp

@@ -96,6 +96,23 @@ ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
    gfi->SyncAliasMemory(*this);
 }
 
+void
+ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff)
+{
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
+   gfr->ProjectCoefficient(real_coeff);
+   *gfi = 0.0;
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
+}
+
+void
+ComplexGridFunction::ProjectCoefficient(ComplexCoefficient &coeff)
+{
+   this->ProjectCoefficient(coeff.real(), coeff.imag());
+}
+
 void
 ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
                                         VectorCoefficient &imag_vcoeff)
@@ -108,6 +125,23 @@ ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
    gfi->SyncAliasMemory(*this);
 }
 
+void
+ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff)
+{
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
+   gfr->ProjectCoefficient(real_vcoeff);
+   *gfi = 0.0;
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
+}
+
+void
+ComplexGridFunction::ProjectCoefficient(ComplexVectorCoefficient &vcoeff)
+{
+   this->ProjectCoefficient(vcoeff.real(), vcoeff.imag());
+}
+
 void
 ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
                                            Coefficient &imag_coeff,
@@ -121,6 +155,26 @@ ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
    gfi->SyncAliasMemory(*this);
 }
 
+void
+ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
+                                           Array<int> &attr)
+{
+   ConstantCoefficient zero_coeff(0.0);
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
+   gfr->ProjectBdrCoefficient(real_coeff, attr);
+   gfi->ProjectBdrCoefficient(zero_coeff, attr);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficient(ComplexCoefficient &coeff,
+                                           Array<int> &attr)
+{
+   this->ProjectBdrCoefficient(coeff.real(), coeff.imag(), attr);
+}
+
 void
 ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
                                                  VectorCoefficient &imag_vcoeff,
@@ -134,6 +188,28 @@ ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
    gfi->SyncAliasMemory(*this);
 }
 
+void
+ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
+                                                 Array<int> &attr)
+{
+   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
+   VectorConstantCoefficient zero_vcoeff(zero_vec);
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
+   gfr->ProjectBdrCoefficientNormal(real_vcoeff, attr);
+   gfi->ProjectBdrCoefficientNormal(zero_vcoeff, attr);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficientNormal(
+   ComplexVectorCoefficient &vcoeff,
+   Array<int> &attr)
+{
+   this->ProjectBdrCoefficientNormal(vcoeff.real(), vcoeff.imag(), attr);
+}
+
 void
 ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
                                                   &real_vcoeff,
@@ -149,6 +225,80 @@ ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
    gfi->SyncAliasMemory(*this);
 }
 
+void
+ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
+                                                  &real_vcoeff,
+                                                  Array<int> &attr)
+{
+   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
+   VectorConstantCoefficient zero_vcoeff(zero_vec);
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
+   gfr->ProjectBdrCoefficientTangent(real_vcoeff, attr);
+   gfi->ProjectBdrCoefficientTangent(zero_vcoeff, attr);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficientTangent(
+   ComplexVectorCoefficient &vcoeff,
+   Array<int> &attr)
+{
+   this->ProjectBdrCoefficientTangent(vcoeff.real(), vcoeff.imag(), attr);
+}
+
+real_t
+ComplexGridFunction::ComputeL2Error(Coefficient &re_exsol,
+                                    Coefficient &im_exsol,
+                                    const IntegrationRule *irs[],
+                                    const Array<int> *elems) const
+{
+   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
+   real_t err_i = gfi->ComputeL2Error(im_exsol, irs, elems);
+
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
+real_t
+ComplexGridFunction::ComputeL2Error(Coefficient &re_exsol,
+                                    const IntegrationRule *irs[],
+                                    const Array<int> *elems) const
+{
+   ConstantCoefficient zero_coef(0.0);
+
+   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
+   real_t err_i = gfi->ComputeL2Error(zero_coef, irs, elems);
+
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
+real_t
+ComplexGridFunction::ComputeL2Error(VectorCoefficient &re_exsol,
+                                    VectorCoefficient &im_exsol,
+                                    const IntegrationRule *irs[],
+                                    const Array<int> *elems) const
+{
+   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
+   real_t err_i = gfi->ComputeL2Error(im_exsol, irs, elems);
+
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
+real_t
+ComplexGridFunction::ComputeL2Error(VectorCoefficient &re_exsol,
+                                    const IntegrationRule *irs[],
+                                    const Array<int> *elems) const
+{
+   Vector zero_vec(re_exsol.GetVDim()); zero_vec = 0.0;
+   VectorConstantCoefficient zero_coef(zero_vec);
+
+   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
+   real_t err_i = gfi->ComputeL2Error(zero_coef, irs, elems);
+
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
 
 ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *fes,
                                      ComplexOperator::Convention convention)
@@ -731,6 +881,17 @@ ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
    pgfi->SyncAliasMemory(*this);
 }
 
+void
+ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff)
+{
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
+   pgfr->ProjectCoefficient(real_coeff);
+   *pgfi = 0.0;
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+}
+
 void
 ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
                                            VectorCoefficient &imag_vcoeff)
@@ -743,6 +904,17 @@ ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
    pgfi->SyncAliasMemory(*this);
 }
 
+void
+ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff)
+{
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
+   pgfr->ProjectCoefficient(real_vcoeff);
+   *pgfi = 0.0;
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+}
+
 void
 ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
                                               Coefficient &imag_coeff,
@@ -756,6 +928,19 @@ ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
    pgfi->SyncAliasMemory(*this);
 }
 
+void
+ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
+                                              Array<int> &attr)
+{
+   ConstantCoefficient zero_coeff(0.0);
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
+   pgfr->ProjectBdrCoefficient(real_coeff, attr);
+   pgfi->ProjectBdrCoefficient(zero_coeff, attr);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+}
+
 void
 ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
                                                     &real_vcoeff,
@@ -771,6 +956,21 @@ ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
    pgfi->SyncAliasMemory(*this);
 }
 
+void
+ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
+                                                    &real_vcoeff,
+                                                    Array<int> &attr)
+{
+   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
+   VectorConstantCoefficient zero_vcoeff(zero_vec);
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
+   pgfr->ProjectBdrCoefficientNormal(real_vcoeff, attr);
+   pgfi->ProjectBdrCoefficientNormal(zero_vcoeff, attr);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+}
+
 void
 ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
                                                      &real_vcoeff,
@@ -786,6 +986,21 @@ ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
    pgfi->SyncAliasMemory(*this);
 }
 
+void
+ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
+                                                     &real_vcoeff,
+                                                     Array<int> &attr)
+{
+   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
+   VectorConstantCoefficient zero_vcoeff(zero_vec);
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
+   pgfr->ProjectBdrCoefficientTangent(real_vcoeff, attr);
+   pgfi->ProjectBdrCoefficientTangent(zero_vcoeff, attr);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+}
+
 void
 ParComplexGridFunction::Distribute(const Vector *tv)
 {
@@ -825,6 +1040,31 @@ ParComplexGridFunction::ParallelProject(Vector &tv) const
    tvi.SyncAliasMemory(tv);
 }
 
+real_t
+ParComplexGridFunction::ComputeL2Error(Coefficient &exsolr,
+                                       const IntegrationRule *irs[],
+                                       Array<int> *elems) const
+{
+   ConstantCoefficient zeroCoef(0.0);
+
+   real_t err_r = pgfr->ComputeL2Error(exsolr, irs, elems);
+   real_t err_i = pgfi->ComputeL2Error(zeroCoef, irs, elems);
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
+real_t
+ParComplexGridFunction::ComputeL2Error(VectorCoefficient &exsolr,
+                                       const IntegrationRule *irs[],
+                                       Array<int> *elems) const
+{
+   Vector zeroVec(exsolr.GetVDim()); zeroVec = 0.0;
+   VectorConstantCoefficient zeroCoef(zeroVec);
+
+   real_t err_r = pgfr->ComputeL2Error(exsolr, irs, elems);
+   real_t err_i = pgfi->ComputeL2Error(zeroCoef, irs, elems);
+   return sqrt(err_r * err_r + err_i * err_i);
+}
+
 
 ParComplexLinearForm::ParComplexLinearForm(ParFiniteElementSpace *pfes,
                                            ComplexOperator::Convention

fem/dfem/doperator.hpp

@@ -241,6 +241,7 @@ public:
    {
       MFEM_ASSERT(!action_callbacks.empty(), "no integrators have been set");
       prolongation(solutions, solutions_t, solutions_l);
+      residual_l = 0.0;
       for (auto &action : action_callbacks)
       {
          action(solutions_l, parameters_l, residual_l);

fem/dfem/util.hpp

@@ -987,7 +987,7 @@ get_restriction_transpose(
    {
       auto RT = [=](const Vector &v_e, Vector &v_l)
       {
-         v_l = v_e;
+         v_l += v_e;
       };
       return std::make_tuple(RT, 1);
    }
@@ -996,7 +996,7 @@ get_restriction_transpose(
       const Operator *R = get_restriction<entity_t>(f, o);
       std::function<void(const Vector&, Vector&)> RT = [=](const Vector &x, Vector &y)
       {
-         R->MultTranspose(x, y);
+         R->AddMultTranspose(x, y);
       };
       return std::make_tuple(RT, R->Height());
    }

general/CMakeLists.txt

@@ -39,6 +39,7 @@ list(APPEND HDRS
   arrays_by_name.hpp
   backends.hpp
   binaryio.hpp
+  complex_type.hpp
   cuda.hpp
   device.hpp
   error.hpp

general/complex_type.hpp

@@ -0,0 +1,125 @@
+// Copyright (c) 2010-2025, 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_COMPLEX_TYPE
+#define MFEM_COMPLEX_TYPE
+
+#include "../config/config.hpp"
+
+#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#include <complex>
+#include <utility>
+#endif
+
+#if defined(MFEM_USE_CUDA)
+#include <cuComplex.h>
+#endif
+
+#if defined(MFEM_USE_HIP)
+#include <hip/hip_complex.h>
+#endif
+
+
+namespace mfem
+{
+
+/// @brief Complex number type for device.
+#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+
+#define zAbs std::abs
+#define zExp std::exp
+#define zNorm std::norm
+using complex_t = std::complex<real_t>;
+
+#else // CUDA or HIP
+
+#if defined(MFEM_USE_CUDA)
+using DoubleComplex_t = cuDoubleComplex;
+#endif
+
+#if defined(MFEM_USE_HIP)
+using DoubleComplex_t = hipDoubleComplex;
+#endif
+
+struct Complex : public DoubleComplex_t
+{
+   MFEM_HOST_DEVICE Complex() = default;
+   MFEM_HOST_DEVICE Complex(real_t r) { x = r, y = 0.0; }
+   MFEM_HOST_DEVICE Complex(real_t r, real_t i) { x = r, y = i; }
+   MFEM_HOST_DEVICE real_t real() const { return x; }
+   MFEM_HOST_DEVICE void real(real_t r) { x = r; }
+   MFEM_HOST_DEVICE real_t imag() const { return y; }
+   MFEM_HOST_DEVICE void imag(real_t i) { y = i; }
+
+   template <typename U>
+   MFEM_HOST_DEVICE inline Complex &operator*=(const U &z)
+   {
+      return *this = *this * z, *this;
+   }
+
+   template <typename U>
+   MFEM_HOST_DEVICE inline Complex &operator/=(const U &z)
+   {
+      return *this = *this / z, *this;
+   }
+};
+
+MFEM_HOST_DEVICE inline Complex operator*(const Complex &x, const real_t &y)
+{
+   return Complex(x.real() * y, x.imag() * y);
+}
+
+MFEM_HOST_DEVICE inline Complex operator+(const Complex &a, const Complex &b)
+{
+   return Complex(a.real() + b.real(), a.imag() + b.imag());
+}
+
+MFEM_HOST_DEVICE inline Complex operator*(const real_t d, const Complex &z)
+{
+   return Complex(z.real() * d, z.imag() * d);
+}
+
+MFEM_HOST_DEVICE inline Complex operator*(const Complex &a, const Complex &b)
+{
+   return Complex(a.real() * b.real() - a.imag() * b.imag(),
+                  a.real() * b.imag() + a.imag() * b.real());
+}
+
+MFEM_HOST_DEVICE inline Complex operator/(const Complex &z, const real_t &d)
+{
+   return Complex(z.real() / d, z.imag() / d);
+}
+
+MFEM_HOST_DEVICE inline real_t zAbs(const Complex &z)
+{
+   return std::hypot(z.real(), z.imag());
+}
+
+MFEM_HOST_DEVICE inline Complex zExp(const Complex &q)
+{
+   Complex z;
+   real_t s, c, e = std::exp(q.real());
+   sincos(q.imag(), &s, &c);
+   z.real(c * e), z.imag(s * e);
+   return z;
+}
+
+MFEM_HOST_DEVICE inline real_t zNorm(const Complex &z)
+{
+   return z.real() * z.real() + z.imag() * z.imag();
+}
+
+using complex_t = Complex;
+#endif // MFEM_USE_CUDA || MFEM_USE_HIP
+
+} // namespace mfem
+
+#endif // MFEM_COMPLEX_TYPE

linalg/CMakeLists.txt

@@ -21,6 +21,7 @@ list(APPEND SRCS
   blockvector.cpp
   complex_densemat.cpp
   complex_operator.cpp
+  complex_vector.cpp
   constraints.cpp
   densemat.cpp
   symmat.cpp
@@ -47,6 +48,7 @@ list(APPEND HDRS
   blockvector.hpp
   complex_densemat.hpp
   complex_operator.hpp
+  complex_vector.hpp
   constraints.hpp
   densemat.hpp
   dinvariants.hpp

linalg/complex_densemat.cpp

@@ -9,6 +9,7 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
+#include "../general/forall.hpp"
 #include "complex_densemat.hpp"
 #include "lapack.hpp"
 #include <complex>
@@ -16,6 +17,8 @@
 namespace mfem
 {
 
+using namespace std;
+
 DenseMatrix & ComplexDenseMatrix::real()
 {
    MFEM_ASSERT(Op_Real_, "ComplexDenseMatrix has no real part!");
@@ -1017,4 +1020,303 @@ void ComplexCholeskyFactors::GetInverseMatrix(int m, real_t * X_r,
    delete [] X;
 }
 
+ComplexTypeDenseMatrix::ComplexTypeDenseMatrix()
+   : height(0), width(0)
+{}
+
+ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(const ComplexTypeDenseMatrix &m)
+   : height(m.Height()), width(m.Width())
+{
+   const int hw = height * width;
+   if (hw > 0)
+   {
+      MFEM_ASSERT(m.data, "invalid source matrix");
+      data.New(hw);
+      std::memcpy(data, m.data, sizeof(complex_t)*hw);
+   }
+}
+
+ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(const DenseMatrix &m)
+   : height(m.Height()), width(m.Width())
+{
+   const int hw = height * width;
+   if (hw > 0)
+   {
+      MFEM_ASSERT(m.data, "invalid source matrix");
+      data.New(hw);
+
+      for (int i = 0; i < hw; i++)
+      {
+         data[i] = m.data[i];
+      }
+   }
+}
+
+ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(int s)
+   : height(s), width(s)
+{
+   MFEM_ASSERT(s >= 0, "invalid DenseMatrix size: " << s);
+   if (s > 0)
+   {
+      data.New(s*s);
+      *this = 0.0; // init with zeroes
+   }
+}
+
+ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(int m, int n)
+   : height(m), width(n)
+{
+   MFEM_ASSERT(m >= 0 && n >= 0,
+               "invalid DenseMatrix size: " << m << " x " << n);
+   const int capacity = m*n;
+   if (capacity > 0)
+   {
+      data.New(capacity);
+      *this = 0.0; // init with zeroes
+   }
+}
+
+void ComplexTypeDenseMatrix::SetSize(int h, int w)
+{
+   MFEM_ASSERT(h >= 0 && w >= 0,
+               "invalid ComplexTypeDenseMatrix size: " << h << " x " << w);
+   if (Height() == h && Width() == w)
+   {
+      return;
+   }
+   height = h;
+   width = w;
+   const int hw = h*w;
+   if (hw > data.Capacity())
+   {
+      data.Delete();
+      data.New(hw);
+      *this = 0.0; // init with zeroes
+   }
+}
+
+/// Returns reference to a_{ij}.
+complex_t &ComplexTypeDenseMatrix::Elem(int i, int j)
+{
+   return (*this)(i,j);
+}
+
+/// Returns constant reference to a_{ij}.
+const complex_t &ComplexTypeDenseMatrix::Elem(int i, int j) const
+{
+   return (*this)(i,j);
+}
+
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(real_t c)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] = c;
+   }
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(complex_t c)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] = c;
+   }
+   return *this;
+}
+
+/// Copy the matrix entries from the given array
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(const real_t *d)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] = d[i];
+   }
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=
+(const complex_t *d)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] = d[i];
+   }
+   return *this;
+}
+
+/// Sets the matrix size and elements equal to those of m
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(const DenseMatrix &m)
+{
+   SetSize(m.height, m.width);
+
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] = m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=
+(const ComplexTypeDenseMatrix &m)
+{
+   SetSize(m.height, m.width);
+
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] = m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=(const real_t *m)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] += m[i];
+   }
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=
+(const complex_t *m)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] += m[i];
+   }
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=(const DenseMatrix &m)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] += m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=
+(const ComplexTypeDenseMatrix &m)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] += m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator-=(const DenseMatrix &m)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] -= m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator-=
+(const ComplexTypeDenseMatrix &m)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] -= m.data[i];
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator*=(real_t c)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] *= c;
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator*=(complex_t c)
+{
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] *= c;
+   }
+
+   return *this;
+}
+
+ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::Set(const DenseMatrix &Mr,
+                                                    const DenseMatrix &Mi)
+{
+   MFEM_ASSERT(height == Mr.Height() && height == Mi.Height() &&
+               width == Mr.Width() && width == Mi.Width(),
+               "incompatible Matrices!");
+
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] = complex_t(Mr.data[i], Mi.data[i]);
+   }
+
+   return *this;
+}
+
+void ComplexTypeDenseMatrix::Swap(ComplexTypeDenseMatrix &other)
+{
+   mfem::Swap(width, other.width);
+   mfem::Swap(height, other.height);
+   mfem::Swap(data, other.data);
+}
+
+ComplexTypeDenseMatrix::~ComplexTypeDenseMatrix()
+{
+   data.Delete();
+}
+
+const DenseMatrix &ComplexTypeDenseMatrix::real() const
+{
+   re_part.SetSize(height, width);
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      re_part.data[i] = data[i].real();
+   }
+
+   return re_part;
+}
+
+const DenseMatrix &ComplexTypeDenseMatrix::imag() const
+{
+   im_part.SetSize(height, width);
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      im_part.data[i] = data[i].imag();
+   }
+
+   return im_part;
+}
+
 } // mfem namespace

linalg/complex_densemat.hpp

@@ -13,6 +13,7 @@
 #define MFEM_COMPLEX_DENSEMAT
 
 #include "complex_operator.hpp"
+#include "../general/complex_type.hpp"
 #include <complex>
 
 namespace mfem
@@ -241,6 +242,220 @@ public:
 
 };
 
+class ComplexTypeDenseMatrix
+{
+protected:
+   int height; ///< Dimension of the output / number of rows in the matrix.
+   int width;  ///< Dimension of the input / number of columns in the matrix.
+
+private:
+   Memory<complex_t > data;
+
+   mutable DenseMatrix re_part;
+   mutable DenseMatrix im_part;
+
+public:
+   /** Default constructor for DenseMatrix.
+       Sets data = NULL and height = width = 0. */
+   ComplexTypeDenseMatrix();
+
+   /// Copy constructor
+   ComplexTypeDenseMatrix(const ComplexTypeDenseMatrix &);
+   ComplexTypeDenseMatrix(const DenseMatrix &);
+
+   /// Creates square matrix of size s.
+   explicit ComplexTypeDenseMatrix(int s);
+
+   /// Creates rectangular matrix of size m x n.
+   ComplexTypeDenseMatrix(int m, int n);
+
+   /// Construct a ComplexTypeDenseMatrix using an existing data array.
+   /** The ComplexTypeDenseMatrix does not assume ownership of the data array,
+       i.e. it will not delete the array. */
+   ComplexTypeDenseMatrix(complex_t *d, int h, int w)
+      : height(h), width(w) { UseExternalData(d, h, w); }
+
+   /// Create a dense matrix using a braced initializer list
+   /// The inner lists correspond to rows of the matrix
+   template <int M, int N, typename T = real_t>
+   explicit ComplexTypeDenseMatrix(const T (&values)[M][N]) :
+      ComplexTypeDenseMatrix(
+         M, N)
+   {
+      // DenseMatrix is column-major so copies have to be element-wise
+      for (int i = 0; i < M; i++)
+      {
+         for (int j = 0; j < N; j++)
+         {
+            (*this)(i,j) = values[i][j];
+         }
+      }
+   }
+
+   /// Change the data array and the size of the DenseMatrix.
+   /** The DenseMatrix does not assume ownership of the data array, i.e. it will
+       not delete the data array @a d. This method should not be used with
+       DenseMatrix that owns its current data array. */
+   void UseExternalData(complex_t *d, int h, int w)
+   {
+      data.Wrap(d, h*w, false);
+      height = h; width = w;
+   }
+
+   /// Change the data array and the size of the DenseMatrix.
+   /** The DenseMatrix does not assume ownership of the data array, i.e. it will
+       not delete the new array @a d. This method will delete the current data
+       array, if owned. */
+   void Reset(complex_t *d, int h, int w)
+   { if (OwnsData()) { data.Delete(); } UseExternalData(d, h, w); }
+
+   /** Clear the data array and the dimensions of the DenseMatrix. This method
+       should not be used with DenseMatrix that owns its current data array. */
+   void ClearExternalData() { data.Reset(); height = width = 0; }
+
+   /// Delete the matrix data array (if owned) and reset the matrix state.
+   void Clear()
+   { if (OwnsData()) { data.Delete(); } ClearExternalData(); }
+
+   /// Get the height (size of output) of the Operator. Synonym with NumRows().
+   inline int Height() const { return height; }
+   /** @brief Get the number of rows (size of output) of the Operator. Synonym
+       with Height(). */
+   inline int NumRows() const { return height; }
+
+   /// Get the width (size of input) of the Operator. Synonym with NumCols().
+   inline int Width() const { return width; }
+   /** @brief Get the number of columns (size of input) of the Operator. Synonym
+       with Width(). */
+   inline int NumCols() const { return width; }
+
+   /// For backward compatibility define Size to be synonym of Width()
+   int Size() const { return Width(); }
+
+   // Total size = width*height
+   int TotalSize() const { return width*height; }
+
+   /// Change the size of the DenseMatrix to s x s.
+   void SetSize(int s) { SetSize(s, s); }
+
+   /// Change the size of the DenseMatrix to h x w.
+   void SetSize(int h, int w);
+
+   /// Returns the matrix data array.
+   inline complex_t *Data() const
+   {
+      return const_cast<complex_t*>
+             ((const complex_t*)data);
+   }
+
+   /// Returns the matrix data array.
+   inline complex_t *GetData() const { return Data(); }
+
+   Memory<complex_t > &GetMemory() { return data; }
+   const Memory<complex_t > &GetMemory() const { return data; }
+
+   /// Return the DenseMatrix data (host pointer) ownership flag.
+   inline bool OwnsData() const { return data.OwnsHostPtr(); }
+
+   /// Returns reference to a_{ij}.
+   inline complex_t &operator()(int i, int j);
+
+   /// Returns constant reference to a_{ij}.
+   inline const complex_t &operator()(int i, int j) const;
+
+   /// Returns reference to a_{ij}.
+   complex_t &Elem(int i, int j);
+
+   /// Returns constant reference to a_{ij}.
+   const complex_t &Elem(int i, int j) const;
+
+   /// Sets the matrix elements equal to constant c
+   ComplexTypeDenseMatrix &operator=(real_t c);
+   ComplexTypeDenseMatrix &operator=(complex_t c);
+
+   /// Copy the matrix entries from the given array
+   ComplexTypeDenseMatrix &operator=(const real_t *d);
+   ComplexTypeDenseMatrix &operator=(const complex_t *d);
+
+   /// Sets the matrix size and elements equal to those of m
+   ComplexTypeDenseMatrix &operator=(const DenseMatrix &m);
+   ComplexTypeDenseMatrix &operator=(const ComplexTypeDenseMatrix &m);
+
+   ComplexTypeDenseMatrix &operator+=(const real_t *m);
+   ComplexTypeDenseMatrix &operator+=(const complex_t *m);
+   ComplexTypeDenseMatrix &operator+=(const DenseMatrix &m);
+   ComplexTypeDenseMatrix &operator+=(const ComplexTypeDenseMatrix &m);
+
+   ComplexTypeDenseMatrix &operator-=(const DenseMatrix &m);
+   ComplexTypeDenseMatrix &operator-=(const ComplexTypeDenseMatrix &m);
+
+   ComplexTypeDenseMatrix &operator*=(real_t c);
+   ComplexTypeDenseMatrix &operator*=(complex_t c);
+
+   /// (*this) = x + i * y
+   ComplexTypeDenseMatrix &Set(const DenseMatrix &x, const DenseMatrix &y);
+
+   std::size_t MemoryUsage() const
+   { return data.Capacity() * sizeof(complex_t); }
+
+   /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
+   const complex_t *Read(bool on_dev = true) const
+   { return mfem::Read(data, Height()*Width(), on_dev); }
+
+   /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
+   const complex_t *HostRead() const
+   { return mfem::Read(data, Height()*Width(), false); }
+
+   /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
+   complex_t *Write(bool on_dev = true)
+   { return mfem::Write(data, Height()*Width(), on_dev); }
+
+   /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
+   complex_t *HostWrite()
+   { return mfem::Write(data, Height()*Width(), false); }
+
+   /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
+   complex_t *ReadWrite(bool on_dev = true)
+   { return mfem::ReadWrite(data, Height()*Width(), on_dev); }
+
+   /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
+   complex_t *HostReadWrite()
+   { return mfem::ReadWrite(data, Height()*Width(), false); }
+
+   void Swap(ComplexTypeDenseMatrix &other);
+
+   /// Return a reference to the real part of this matrix
+   const DenseMatrix &real() const;
+
+   /// Return a reference to the imaginary part of this matrix
+   const DenseMatrix &imag() const;
+
+   /// Destroys dense matrix.
+   virtual ~ComplexTypeDenseMatrix();
+};
+
+/// Specialization of the template function Swap<> for class ComplexTypeDenseMatrix
+template<> inline void Swap<ComplexTypeDenseMatrix>(ComplexTypeDenseMatrix &a,
+                                                    ComplexTypeDenseMatrix &b)
+{
+   a.Swap(b);
+}
+
+// Inline methods
+
+inline complex_t &ComplexTypeDenseMatrix::operator()(int i, int j)
+{
+   MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
+   return data[i+j*height];
+}
+
+inline const complex_t &ComplexTypeDenseMatrix::operator()
+(int i, int j) const
+{
+   MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
+   return data[i+j*height];
+}
+
 } // namespace mfem
 
 #endif // MFEM_COMPLEX_DENSEMAT

linalg/complex_vector.cpp

@@ -0,0 +1,424 @@
+// Copyright (c) 2010-2025, 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 "../general/forall.hpp"
+#include "../general/reducers.hpp"
+#include "complex_vector.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+ComplexVector::ComplexVector(const ComplexVector &v)
+{
+   const int s = v.Size();
+   size = s;
+   if (s > 0)
+   {
+      MFEM_ASSERT(!v.data.Empty(), "invalid source vector");
+      data.New(s, v.data.GetMemoryType());
+      data.CopyFrom(v.data, s);
+   }
+   UseDevice(v.UseDevice());
+}
+
+ComplexVector::ComplexVector(const Vector &v)
+{
+   const int s = v.Size();
+   size = s;
+   if (s > 0)
+   {
+      MFEM_ASSERT(!v.data.Empty(), "invalid source vector");
+      data.New(s, v.data.GetMemoryType());
+      MFEM_FORALL(i, size, data[i] = v.data[i]; );
+   }
+   UseDevice(v.UseDevice());
+}
+
+ComplexVector::ComplexVector(ComplexVector &&v)
+{
+   *this = std::move(v);
+}
+
+complex_t &ComplexVector::Elem(int i)
+{
+   return operator()(i);
+}
+
+const complex_t &ComplexVector::Elem(int i) const
+{
+   return operator()(i);
+}
+
+complex_t ComplexVector::operator*(const complex_t *v) const
+{
+   HostRead();
+   complex_t dot = 0.0;
+#ifdef MFEM_USE_LEGACY_OPENMP
+   #pragma omp parallel for reduction(+:dot)
+#endif
+   for (int i = 0; i < size; i++)
+   {
+      dot += data[i] * v[i];
+   }
+   return dot;
+}
+
+complex_t ComplexVector::operator*(const real_t *v) const
+{
+   HostRead();
+   complex_t dot = 0.0;
+#ifdef MFEM_USE_LEGACY_OPENMP
+   #pragma omp parallel for reduction(+:dot)
+#endif
+   for (int i = 0; i < size; i++)
+   {
+      dot += data[i] * v[i];
+   }
+   return dot;
+}
+
+complex_t ComplexVector::operator*(const ComplexVector &v) const
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   if (size == 0) { return 0.0; }
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const auto m_data = Read(use_dev), v_data = v.Read(use_dev);
+
+   // The standard way of computing the dot product is non-deterministic
+   complex_t prod = 0.0;
+   for (int i = 0; i < size; i++)
+   {
+      prod += m_data[i] * v_data[i];
+   }
+   return prod;
+}
+
+complex_t ComplexVector::operator*(const Vector &v) const
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   if (size == 0) { return 0.0; }
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const auto m_data = Read(use_dev);
+   const auto v_data = v.Read(use_dev);
+
+   // The standard way of computing the dot product is non-deterministic
+   complex_t prod = 0.0;
+   for (int i = 0; i < size; i++)
+   {
+      prod += m_data[i] * v_data[i];
+   }
+   return prod;
+}
+
+ComplexVector &ComplexVector::operator=(const complex_t *v)
+{
+   HostRead();
+   MFEM_FORALL(i, size, data[i] = v[i]; );
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(const real_t *v)
+{
+   HostRead();
+   MFEM_FORALL(i, size, data[i] = v[i]; );
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(const ComplexVector &v)
+{
+#if 0
+   SetSize(v.Size(), v.data.GetMemoryType());
+   data.CopyFrom(v.data, v.Size());
+   UseDevice(v.UseDevice());
+#else
+   SetSize(v.Size());
+   const bool vuse = v.UseDevice();
+   const bool use_dev = UseDevice() || vuse;
+   v.UseDevice(use_dev);
+   // keep 'data' where it is, unless 'use_dev' is true
+   if (use_dev) { Write(); }
+   data.CopyFrom(v.data, v.Size());
+   v.UseDevice(vuse);
+#endif
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(const Vector &v)
+{
+   SetSize(v.Size());
+   const bool vuse = v.UseDevice();
+   const bool use_dev = UseDevice() || vuse;
+   v.UseDevice(use_dev);
+   // keep 'data' where it is, unless 'use_dev' is true
+   if (use_dev) { Write(); }
+   MFEM_FORALL(i, size, data[i] = v[i]; );
+   v.UseDevice(vuse);
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(ComplexVector &&v)
+{
+   v.Swap(*this);
+   if (this != &v) { v.Destroy(); }
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(complex_t value)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = Write(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] = value; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator=(real_t value)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = Write(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] = value; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator*=(complex_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator*=(real_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator*=(const ComplexVector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator*=(const Vector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator/=(complex_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   const complex_t m = conj(c) / norm(c);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= m; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator/=(real_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   const real_t m = 1.0/c;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] *= m; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator/=(const ComplexVector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] /= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator/=(const Vector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] /= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator-=(complex_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] -= c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator-=(real_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] -= c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator-=(const ComplexVector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] -= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator-=(const Vector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] -= x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator+=(complex_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] += c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator+=(real_t c)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] += c; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator+=(const ComplexVector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] += x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::operator+=(const Vector &v)
+{
+   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || v.UseDevice();
+   const int N = size;
+   const auto x = v.Read(use_dev);
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] += x[i]; });
+   return *this;
+}
+
+ComplexVector &ComplexVector::Set(const Vector &Vr, const Vector &Vi)
+{
+   MFEM_ASSERT(size == Vr.size && size == Vi.size, "incompatible Vectors!");
+
+   const bool use_dev = UseDevice() || Vr.UseDevice() || Vi.UseDevice();
+   const int N = size;
+   const auto x = Vr.Read(use_dev);
+   const auto y = Vi.Read(use_dev);
+   auto z = Write(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { z[i] = complex_t(x[i], y[i]); });
+   return *this;
+}
+
+const Vector &ComplexVector::real() const
+{
+   re_part.SetSize(size);
+   const bool use_dev = UseDevice();
+   const int N = size;
+   const auto z = Read(use_dev);
+   auto x = re_part.Write(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { x[i] = z[i].real(); });
+   return re_part;
+}
+
+const Vector &ComplexVector::imag() const
+{
+   im_part.SetSize(size);
+   const bool use_dev = UseDevice();
+   const int N = size;
+   const auto z = Read(use_dev);
+   auto y = im_part.Write(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   { y[i] = z[i].imag(); });
+   return im_part;
+}
+
+}

linalg/complex_vector.hpp

@@ -0,0 +1,479 @@
+// Copyright (c) 2010-2025, 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_COMPLEX_VECTOR
+#define MFEM_COMPLEX_VECTOR
+
+#include "vector.hpp"
+#include "../general/complex_type.hpp"
+
+namespace mfem
+{
+
+class ComplexVector
+{
+private:
+
+   Memory<complex_t > data;
+   int size;
+
+   mutable Vector re_part;
+   mutable Vector im_part;
+
+public:
+
+   /// Default constructor for ComplexVector. Sets size = 0
+   ComplexVector() : size(0) { }
+
+   /// Copy constructor. Allocates a new data array and copies the data.
+   ComplexVector(const ComplexVector &);
+
+   /// Copy constructor. Allocates a new data array and copies the
+   /// data into real part of this vector.
+   ComplexVector(const Vector &);
+
+   /// Move constructor. "Steals" data from its argument.
+   ComplexVector(ComplexVector&& v);
+
+   /// @brief Creates vector of size s.
+   /// @warning Entries are not initialized to zero!
+   explicit ComplexVector(int s);
+
+   /// Creates a vector referencing an array of complex<doubles>,
+   /// owned by someone else.
+   /// The pointer @a data_ can be NULL. The data array can be replaced later
+   /// with SetData().
+   ComplexVector(complex_t *data_, int size_)
+   { data.Wrap(data_, size_, false); size = size_; }
+
+   /// @brief Create a ComplexVector referencing a sub-vector of the
+   //  ComplexVector @a base starting at the given offset, @a
+   //  base_offset, and size @a size_.
+   ComplexVector(ComplexVector &base, int base_offset, int size_)
+      : data(base.data, base_offset, size_), size(size_) { }
+
+   /// Create a ComplexVector of size @a size_ using MemoryType @a mt.
+   ComplexVector(int size_, MemoryType mt)
+      : data(size_, mt), size(size_) { }
+
+   /// @brief Create a ComplexVector of size @a size_ using host
+   /// MemoryType @a h_mt and device MemoryType @a d_mt.
+   ComplexVector(int size_, MemoryType h_mt, MemoryType d_mt)
+      : data(size_, h_mt, d_mt), size(size_) { }
+
+   /// Create a vector from a statically sized C-style array of convertible type
+   template <typename CT, int N>
+   explicit ComplexVector(const CT (&values)[N]) : ComplexVector(N)
+   { std::copy(values, values + N, begin()); }
+
+   /// Create a vector using a braced initializer list
+   template <typename CT, typename std::enable_if<
+                std::is_convertible<CT,complex_t >::value,bool>::type = true>
+   explicit ComplexVector(std::initializer_list<CT> values) : ComplexVector(
+         values.size())
+   { std::copy(values.begin(), values.end(), begin()); }
+
+   /// Enable execution of Vector operations using the mfem::Device.
+   /// The default is to use Backend::CPU (serial execution on each MPI rank),
+   /// regardless of the mfem::Device configuration.
+   ///
+   /// When appropriate, MFEM functions and class methods will enable the use
+   /// of the mfem::Device for their Vector parameters.
+   ///
+   /// Some derived classes, e.g. GridFunction, enable the use of the
+   /// mfem::Device by default.
+   virtual void UseDevice(bool use_dev) const { data.UseDevice(use_dev); }
+
+   /// Return the device flag of the Memory object used by the Vector
+   virtual bool UseDevice() const { return data.UseDevice(); }
+
+   /// @brief Resize the vector to size @a s.
+   /// If the new size is less than or equal to Capacity() then the internal
+   /// data array remains the same. Otherwise, the old array is deleted, if
+   /// owned, and a new array of size @a s is allocated without copying the
+   /// previous content of the ComplexVector.
+   /// @warning In the second case above (new size greater than current one),
+   /// the vector will allocate new data array, even if it did not own the
+   /// original data! Also, new entries are not initialized!
+   void SetSize(int s);
+
+   /// Resize the vector to size @a s using MemoryType @a mt.
+   void SetSize(int s, MemoryType mt);
+
+   /// Resize the vector to size @a s using the MemoryType of @a v.
+   void SetSize(int s, const ComplexVector &v)
+   { SetSize(s, v.GetMemory().GetMemoryType()); }
+
+   /// Resize the vector to size @a s using the MemoryType of @a v.
+   void SetSize(int s, const Vector &v)
+   { SetSize(s, v.GetMemory().GetMemoryType()); }
+
+   /// Set the Vector data.
+   /// @warning This method should be called only when OwnsData() is false.
+   void SetData(complex_t *d)
+   { data.Wrap(d, data.Capacity(), false); }
+
+   /// Set the Vector data and size.
+   /// The Vector does not assume ownership of the new data. The new size is
+   /// also used as the new Capacity().
+   /// @warning This method should be called only when OwnsData() is false.
+   /// @sa NewDataAndSize().
+   void SetDataAndSize(complex_t *d, int s)
+   { data.Wrap(d, s, false); size = s; }
+
+   /// Set the Vector data and size, deleting the old data, if owned.
+   /// The Vector does not assume ownership of the new data. The new size is
+   /// also used as the new Capacity().
+   /// @sa SetDataAndSize().
+   void NewDataAndSize(complex_t *d, int s)
+   {
+      data.Delete();
+      SetDataAndSize(d, s);
+   }
+
+   /// Reset the Vector to use the given external Memory @a mem and size @a s.
+   /// If @a own_mem is false, the Vector will not own any of the pointers of
+   /// @a mem.
+   ///
+   /// Note that when @a own_mem is true, the @a mem object can be destroyed
+   /// immediately by the caller but `mem.Delete()` should NOT be called since
+   /// the Vector object takes ownership of all pointers owned by @a mem.
+   ///
+   /// @sa NewDataAndSize().
+   inline void NewMemoryAndSize(const Memory<complex_t > &mem,
+                                int s, bool own_mem);
+
+   /// Reset the Vector to be a reference to a sub-vector of @a base.
+   inline void MakeRef(ComplexVector &base, int offset, int size);
+
+   /// @brief Reset the Vector to be a reference to a sub-vector of @a base
+   /// without changing its current size.
+   inline void MakeRef(ComplexVector &base, int offset);
+
+   /// Set the Vector data (host pointer) ownership flag.
+   void MakeDataOwner() const { data.SetHostPtrOwner(true); }
+
+   /// Destroy a vector
+   void Destroy();
+
+   /// @brief Delete the device pointer, if owned. If @a copy_to_host is true
+   /// and the data is valid only on device, move it to host before deleting.
+   /// Invalidates the device memory.
+   void DeleteDevice(bool copy_to_host = true)
+   { data.DeleteDevice(copy_to_host); }
+
+   /// Returns the size of the vector.
+   inline int Size() const { return size; }
+
+   /// Return the size of the currently allocated data array.
+   /// It is always true that Capacity() >= Size().
+   inline int Capacity() const { return data.Capacity(); }
+
+   /// Return a pointer to the beginning of the ComplexVector data.
+   /// @warning This method should be used with caution as it gives write access
+   /// to the data of const-qualified ComplexVector%s.
+   inline complex_t *GetData() const
+   { return const_cast<complex_t*>((const complex_t*)data); }
+
+   /// STL-like begin.
+   inline complex_t *begin() { return data; }
+
+   /// STL-like end.
+   inline complex_t *end() { return data + size; }
+
+   /// STL-like begin (const version).
+   inline const complex_t *begin() const { return data; }
+
+   /// STL-like end (const version).
+   inline const complex_t *end() const { return data + size; }
+
+   /// Return a reference to the Memory object used by the Vector.
+   Memory<complex_t > &GetMemory() { return data; }
+
+   /// @brief Return a reference to the Memory object used by the
+   /// ComplexVector, const version.
+   const Memory<complex_t > &GetMemory() const { return data; }
+
+   /// Update the memory location of the vector to match @a v.
+   void SyncMemory(const ComplexVector &v) const
+   { GetMemory().Sync(v.GetMemory()); }
+
+   /// Update the alias memory location of the vector to match @a v.
+   void SyncAliasMemory(const ComplexVector &v) const
+   { GetMemory().SyncAlias(v.GetMemory(),Size()); }
+
+   /// Read the Vector data (host pointer) ownership flag.
+   inline bool OwnsData() const { return data.OwnsHostPtr(); }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline void StealData(complex_t **p)
+   { *p = data; data.Reset(); size = 0; }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline complex_t *StealData()
+   { complex_t *p; StealData(&p); return p; }
+
+   /// Access Vector entries. Index i = 0 .. size-1.
+   complex_t &Elem(int i);
+
+   /// Read only access to Vector entries. Index i = 0 .. size-1.
+   const complex_t &Elem(int i) const;
+
+   /// Access Vector entries using () for 0-based indexing.
+   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
+   inline complex_t &operator()(int i);
+
+   /// Read only access to Vector entries using () for 0-based indexing.
+   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
+   inline const complex_t &operator()(int i) const;
+
+   /// Access Vector entries using [] for 0-based indexing.
+   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
+   inline complex_t &operator[](int i) { return (*this)(i); }
+
+   /// Read only access to Vector entries using [] for 0-based indexing.
+   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
+   inline const complex_t &operator[](int i) const
+   { return (*this)(i); }
+
+   /// Dot product with a `complex<double> *` array.
+   /// @note No complex conjugate is performed
+   complex_t operator*(const complex_t *v) const;
+   complex_t operator*(const real_t *v) const;
+
+   /// Return the inner-product.
+   /// @note No complex conjugate is performed
+   complex_t operator*(const ComplexVector &v) const;
+   complex_t operator*(const Vector &v) const;
+
+   /// Copy Size() entries from @a v.
+   ComplexVector &operator=(const complex_t *v);
+   ComplexVector &operator=(const real_t *v);
+
+   /// Copy assignment.
+   /// @note Defining this method overwrites the implicitly defined copy
+   /// assignment operator.
+   ComplexVector &operator=(const ComplexVector &v);
+   ComplexVector &operator=(const Vector &v);
+
+   /// Move assignment
+   ComplexVector &operator=(ComplexVector&& v);
+
+   /// Redefine '=' for vector = constant.
+   ComplexVector &operator=(complex_t value);
+   ComplexVector &operator=(real_t value);
+
+   /// Scale vector by a constant
+   ComplexVector &operator*=(complex_t c);
+   ComplexVector &operator*=(real_t c);
+
+   /// Component-wise scaling: (*this)(i) *= v(i)
+   ComplexVector &operator*=(const ComplexVector &v);
+   ComplexVector &operator*=(const Vector &v);
+
+   /// Divide vector by a consant
+   ComplexVector &operator/=(complex_t c);
+   ComplexVector &operator/=(real_t c);
+
+   /// Component-wise division: (*this)(i) /= v(i)
+   ComplexVector &operator/=(const ComplexVector &v);
+   ComplexVector &operator/=(const Vector &v);
+
+   /// Subtract a constant from this vector
+   ComplexVector &operator-=(complex_t c);
+   ComplexVector &operator-=(real_t c);
+
+   /// Subtract a vector from this vector
+   ComplexVector &operator-=(const ComplexVector &v);
+   ComplexVector &operator-=(const Vector &v);
+
+   /// Add a constant to this vector
+   ComplexVector &operator+=(complex_t c);
+   ComplexVector &operator+=(real_t c);
+
+   /// Add a vector to this vector
+   ComplexVector &operator+=(const ComplexVector &v);
+   ComplexVector &operator+=(const Vector &v);
+
+   /// (*this) = x + i * y
+   ComplexVector &Set(const Vector &x, const Vector &y);
+
+   /// Swap the contents of two Vectors
+   inline void Swap(ComplexVector &other);
+
+   /// Return a reference to the real part of this vector
+   const Vector &real() const;
+
+   /// Return a reference to the imaginary part of this vector
+   const Vector &imag() const;
+
+   /// Destroys vector.
+   virtual ~ComplexVector();
+
+   /// Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
+   virtual const complex_t *Read(bool on_dev = true) const
+   { return mfem::Read(data, size, on_dev); }
+
+   /// Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
+   virtual const complex_t *HostRead() const
+   { return mfem::Read(data, size, false); }
+
+   /// Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
+   virtual complex_t *Write(bool on_dev = true)
+   { return mfem::Write(data, size, on_dev); }
+
+   /// Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
+   virtual complex_t *HostWrite()
+   { return mfem::Write(data, size, false); }
+
+   /// Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
+   virtual complex_t *ReadWrite(bool on_dev = true)
+   { return mfem::ReadWrite(data, size, on_dev); }
+
+   /// Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
+   virtual complex_t *HostReadWrite()
+   { return mfem::ReadWrite(data, size, false); }
+};
+
+inline ComplexVector::ComplexVector(int s)
+{
+   MFEM_ASSERT(s>=0,"Unexpected negative size.");
+   size = s;
+   if (s > 0)
+   {
+      data.New(s);
+   }
+}
+
+inline void ComplexVector::SetSize(int s)
+{
+   if (s == size)
+   {
+      return;
+   }
+   if (s <= data.Capacity())
+   {
+      size = s;
+      return;
+   }
+   // preserve a valid MemoryType and device flag
+   const MemoryType mt = data.GetMemoryType();
+   const bool use_dev = data.UseDevice();
+   data.Delete();
+   size = s;
+   data.New(s, mt);
+   data.UseDevice(use_dev);
+}
+
+inline void ComplexVector::SetSize(int s, MemoryType mt)
+{
+   if (mt == data.GetMemoryType())
+   {
+      if (s == size)
+      {
+         return;
+      }
+      if (s <= data.Capacity())
+      {
+         size = s;
+         return;
+      }
+   }
+   const bool use_dev = data.UseDevice();
+   data.Delete();
+   if (s > 0)
+   {
+      data.New(s, mt);
+      size = s;
+   }
+   else
+   {
+      data.Reset();
+      size = 0;
+   }
+   data.UseDevice(use_dev);
+}
+
+inline void ComplexVector::NewMemoryAndSize(
+   const Memory<complex_t > &mem,
+   int s,
+   bool own_mem)
+{
+   data.Delete();
+   size = s;
+   if (own_mem)
+   {
+      data = mem;
+   }
+   else
+   {
+      data.MakeAlias(mem, 0, s);
+   }
+}
+
+inline void ComplexVector::MakeRef(ComplexVector &base, int offset, int s)
+{
+   data.Delete();
+   size = s;
+   data.MakeAlias(base.GetMemory(), offset, s);
+}
+
+inline void ComplexVector::MakeRef(ComplexVector &base, int offset)
+{
+   data.Delete();
+   data.MakeAlias(base.GetMemory(), offset, size);
+}
+
+inline void ComplexVector::Destroy()
+{
+   const bool use_dev = data.UseDevice();
+   data.Delete();
+   size = 0;
+   data.Reset();
+   data.UseDevice(use_dev);
+}
+
+inline complex_t &ComplexVector::operator()(int i)
+{
+   MFEM_ASSERT(data && i >= 0 && i < size,
+               "index [" << i << "] is out of range [0," << size << ")");
+
+   return data[i];
+}
+
+inline const complex_t &ComplexVector::operator()(int i) const
+{
+   MFEM_ASSERT(data && i >= 0 && i < size,
+               "index [" << i << "] is out of range [0," << size << ")");
+
+   return data[i];
+}
+
+inline void ComplexVector::Swap(ComplexVector &other)
+{
+   mfem::Swap(data, other.data);
+   mfem::Swap(size, other.size);
+}
+
+/// Specialization of the template function Swap<> for class ComplexVector
+template<> inline void Swap<ComplexVector>(ComplexVector &a, ComplexVector &b)
+{
+   a.Swap(b);
+}
+
+inline ComplexVector::~ComplexVector()
+{
+   data.Delete();
+}
+
+} // namespace mfem
+
+#endif

linalg/densemat.cpp

@@ -4405,4 +4405,32 @@ void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X)
    BatchedLinAlg::LUSolve(Mlu, P, X);
 }
 
+#ifdef MFEM_USE_LAPACK
+void BandedSolve(int KL, int KU, DenseMatrix &AB, DenseMatrix &B,
+                 Array<int> &ipiv)
+{
+   int LDAB = (2*KL) + KU + 1;
+   int N = AB.NumCols();
+   int NRHS = B.NumCols();
+   int info;
+   ipiv.SetSize(N);
+   MFEM_LAPACK_PREFIX(gbsv_)(&N, &KL, &KU, &NRHS, AB.GetData(), &LDAB,
+                             ipiv.GetData(), B.GetData(), &N, &info);
+   MFEM_ASSERT(info == 0, "BandedSolve failed in LAPACK");
+}
+
+void BandedFactorizedSolve(int KL, int KU, DenseMatrix &AB, DenseMatrix &B,
+                           bool transpose, Array<int> &ipiv)
+{
+   int LDAB = (2*KL) + KU + 1;
+   int N = AB.NumCols();
+   int NRHS = B.NumCols();
+   char trans = transpose ? 'T' : 'N';
+   int info;
+   MFEM_LAPACK_PREFIX(gbtrs_)(&trans, &N, &KL, &KU, &NRHS, AB.GetData(), &LDAB,
+                              ipiv.GetData(), B.GetData(), &N, &info);
+   MFEM_ASSERT(info == 0, "BandedFactorizedSolve failed in LAPACK");
+}
+#endif
+
 } // namespace mfem

linalg/densemat.hpp

@@ -24,6 +24,7 @@ class DenseMatrix : public Matrix
 {
    friend class DenseTensor;
    friend class DenseMatrixInverse;
+   friend class ComplexTypeDenseMatrix;
 
 private:
    Memory<real_t> data;
@@ -1329,6 +1330,13 @@ void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const real_t TOL = 0.0);
     dimension m x n. */
 void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X);
 
+#ifdef MFEM_USE_LAPACK
+void BandedSolve(int KL, int KU, DenseMatrix &AB, DenseMatrix &B,
+                 Array<int> &ipiv);
+void BandedFactorizedSolve(int KL, int KU, DenseMatrix &AB, DenseMatrix &B,
+                           bool transpose, Array<int> &ipiv);
+#endif
+
 // Inline methods
 
 inline real_t &DenseMatrix::operator()(int i, int j)

linalg/lapack.hpp

@@ -42,6 +42,13 @@ extern "C" void
 MFEM_LAPACK_PREFIX(getri_)(int *N, real_t *A, int *LDA, int *IPIV, real_t *WORK,
                            int *LWORK, int *INFO);
 extern "C" void
+MFEM_LAPACK_PREFIX(gbsv_)(int *, int *, int *, int *, real_t *, int *, int *,
+                          real_t *, int *, int *);
+extern "C" void
+MFEM_LAPACK_PREFIX(gbtrs_)(char *, int *, int *, int *, int *, real_t *, int *,
+                           int *, real_t *, int *, int *);
+
+extern "C" void
 MFEM_LAPACK_PREFIX(syevr_)(char *JOBZ, char *RANGE, char *UPLO, int *N,
                            real_t *A, int *LDA, real_t *VL, real_t *VU, int *IL,
                            int *IU, real_t *ABSTOL, int *M, real_t *W,

linalg/vector.hpp

@@ -80,6 +80,8 @@ inline real_t rand_real()
 /// Vector data type.
 class Vector
 {
+   friend class ComplexVector;
+
 protected:
 
    Memory<real_t> data;

mesh/mesh.hpp

@@ -588,9 +588,10 @@ protected:
    void Loader(std::istream &input, int generate_edges = 0,
                std::string parse_tag = "");
 
-   /** If NURBS mesh, write NURBS format. If NCMesh, write mfem v1.1 format.
-       If section_delimiter is empty, write mfem v1.0 format. Otherwise, write
-       mfem v1.2 format with the given section_delimiter at the end.
+   /** @brief If NURBS mesh, write NURBS format. If NCMesh, write mfem v1.1
+       format. If section_delimiter is empty, write mfem v1.0 format. Otherwise,
+       write mfem v1.2 format with the given section_delimiter at the end.
+
        If @a comments is non-empty, it will be printed after the first line of
        the file, and each line should begin with '#'. */
    void Printer(std::ostream &os = mfem::out,
@@ -2482,10 +2483,12 @@ public:
    /// Print the mesh to the given stream using Netgen/Truegrid format.
    virtual void PrintXG(std::ostream &os = mfem::out) const;
 
-   /// Print the mesh to the given stream using the default MFEM mesh format.
-   /// \see mfem::ofgzstream() for on-the-fly compression of ascii outputs. If
-   /// @a comments is non-empty, it will be printed after the first line of the
-   /// file, and each line should begin with '#'.
+   /** @brief Print the mesh to the given stream using the default MFEM mesh
+       format.
+
+       \see mfem::ofgzstream() for on-the-fly compression of ascii outputs. If
+       @a comments is non-empty, it will be printed after the first line of the
+       file, and each line should begin with '#'. */
    virtual void Print(std::ostream &os = mfem::out,
                       const std::string &comments = "") const
    { Printer(os, "", comments); }

mesh/nurbs.cpp

@@ -53,7 +53,7 @@ KnotVector::KnotVector(int order, int NCP)
 }
 
 KnotVector::KnotVector(int order, const Vector& intervals,
-                       const Array<int>& continuity )
+                       const Array<int>& continuity)
 {
    // NOTE: This may need to be generalized to support periodicity
    // in the future.
@@ -151,7 +151,7 @@ void KnotVector::UniformRefinement(Vector &newknots, int rf) const
       {
          for (int m = 1; m < rf; ++m)
          {
-            newknots(j) = m * h * (knot(i) + knot(i+1));
+            newknots(j) = ((1.0 - (m * h)) * knot(i)) + (m * h * knot(i+1));
             j++;
          }
       }
@@ -340,7 +340,7 @@ void KnotVector::PrintFunctions(std::ostream &os, int samples) const
    }
 }
 
-// Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
+// Routine from "The NURBS Book" - 2nd ed - Piegl and Tiller
 // Algorithm A2.2 p. 70
 void KnotVector::CalcShape(Vector &shape, int i, real_t xi) const
 {
@@ -367,7 +367,7 @@ void KnotVector::CalcShape(Vector &shape, int i, real_t xi) const
    }
 }
 
-// Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
+// Routine from "The NURBS Book" - 2nd ed - Piegl and Tiller
 // Algorithm A2.3 p. 72
 void KnotVector::CalcDShape(Vector &grad, int i, real_t xi) const
 {
@@ -425,7 +425,7 @@ void KnotVector::CalcDShape(Vector &grad, int i, real_t xi) const
    }
 }
 
-// Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
+// Routine from "The NURBS Book" - 2nd ed - Piegl and Tiller
 // Algorithm A2.3 p. 72
 void KnotVector::CalcDnShape(Vector &gradn, int n, int i, real_t xi) const
 {
@@ -545,11 +545,11 @@ void KnotVector::FindMaxima(Array<int> &ks, Vector &xi, Vector &u) const
          int i = j - d;
          if (isElement(i))
          {
-            arg1 = 1e-16;
+            arg1 = std::numeric_limits<real_t>::epsilon() / 2_r;
             CalcShape(shape, i, arg1);
             max1 = shape[d];
 
-            arg2 = 1-(1e-16);
+            arg2 = 1_r - arg1;
             CalcShape(shape, i, arg2);
             max2 = shape[d];
 
@@ -587,9 +587,9 @@ void KnotVector::FindMaxima(Array<int> &ks, Vector &xi, Vector &u) const
    }
 }
 
-// Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
+// Routine from "The NURBS Book" - 2nd ed - Piegl and Tiller
 // Algorithm A9.1 p. 369
-void KnotVector::FindInterpolant(Array<Vector*> &x)
+void KnotVector::FindInterpolant(Array<Vector*> &x, bool reuse_inverse)
 {
    int order = GetOrder();
    int ncp = GetNCP();
@@ -597,29 +597,93 @@ void KnotVector::FindInterpolant(Array<Vector*> &x)
    // Find interpolation points
    Vector xi_args, u_args;
    Array<int> i_args;
-   FindMaxima(i_args,xi_args, u_args);
+   FindMaxima(i_args, xi_args, u_args);
 
    // Assemble collocation matrix
-   Vector shape(order+1);
-   DenseMatrix A(ncp,ncp);
-   A = 0.0;
+#ifdef MFEM_USE_LAPACK
+   // If using LAPACK, we use banded matrix storage (order + 1 nonzeros per row).
+   // Find banded structure of matrix.
+   int KL = 0; // Number of subdiagonals
+   int KU = 0; // Number of superdiagonals
    for (int i = 0; i < ncp; i++)
    {
-      CalcShape(shape, i_args[i], xi_args[i]);
       for (int p = 0; p < order+1; p++)
       {
-         A(i,i_args[i] + p) =  shape[p];
+         const int col = i_args[i] + p;
+         if (col < i)
+         {
+            KL = std::max(KL, i - col);
+         }
+         else if (i < col)
+         {
+            KU = std::max(KU, col - i);
+         }
       }
    }
 
-   // Solve problems
-   A.Invert();
+   const int LDAB = (2*KL) + KU + 1;
+   const int N = ncp;
+
+   fact_AB.SetSize(LDAB, N);
+#else
+   // Without LAPACK, we store and invert a DenseMatrix (inefficient).
+   if (!reuse_inverse)
+   {
+      A_coll_inv.SetSize(ncp, ncp);
+      A_coll_inv = 0.0;
+   }
+#endif
+
+   Vector shape(order+1);
+
+   if (!reuse_inverse) // Set collocation matrix entries
+   {
+      for (int i = 0; i < ncp; i++)
+      {
+         CalcShape(shape, i_args[i], xi_args[i]);
+         for (int p = 0; p < order+1; p++)
+         {
+            const int j = i_args[i] + p;
+#ifdef MFEM_USE_LAPACK
+            fact_AB(KL+KU+i-j,j) = shape[p];
+#else
+            A_coll_inv(i,j) = shape[p];
+#endif
+         }
+      }
+   }
+
+   // Solve the system
+#ifdef MFEM_USE_LAPACK
+   const int NRHS = x.Size();
+   DenseMatrix B(N, NRHS);
+   for (int j=0; j<NRHS; ++j)
+   {
+      for (int i=0; i<N; ++i) { B(i, j) = (*x[j])[i]; }
+   }
+
+   if (reuse_inverse)
+   {
+      BandedFactorizedSolve(KL, KU, fact_AB, B, false, fact_ipiv);
+   }
+   else
+   {
+      BandedSolve(KL, KU, fact_AB, B, fact_ipiv);
+   }
+
+   for (int j=0; j<NRHS; ++j)
+   {
+      for (int i=0; i<N; ++i) { (*x[j])[i] = B(i, j); }
+   }
+#else
+   if (!reuse_inverse) { A_coll_inv.Invert(); }
    Vector tmp;
-   for (int i= 0; i < x.Size(); i++)
+   for (int i = 0; i < x.Size(); i++)
    {
       tmp = *x[i];
-      A.Mult(tmp,*x[i]);
+      A_coll_inv.Mult(tmp, *x[i]);
    }
+#endif
 }
 
 int KnotVector::findKnotSpan(real_t u) const
@@ -1421,7 +1485,7 @@ void NURBSPatch::DegreeElevate(int t)
    }
 }
 
-// Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
+// Routine from "The NURBS Book" - 2nd ed - Piegl and Tiller
 void NURBSPatch::DegreeElevate(int dir, int t)
 {
    if (dir >= kv.Size() || dir < 0)
@@ -2385,7 +2449,7 @@ NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
 }
 
 NURBSExtension::NURBSExtension(const Mesh *patch_topology,
-                               const Array<const NURBSPatch*> patches_)
+                               const Array<const NURBSPatch*> &patches_)
 {
    // Basic topology checks
    MFEM_VERIFY(patches_.Size() > 0, "Must have at least one patch");

mesh/nurbs.hpp

@@ -143,8 +143,13 @@ public:
    /** @brief Global curve interpolation through the points @a x (overwritten).
        @a x is an array with the length of the spatial dimension containing
        vectors with spatial coordinates. The control points of the interpolated
-       curve are returned in @a x in the same form. */
-   void FindInterpolant(Array<Vector*> &x);
+       curve are returned in @a x in the same form.
+
+       The inverse of the collocation matrix, used in the interpolation, is
+       stored for repeated calls and used if @a reuse_inverse is true. Reuse is
+       valid only if this KnotVector has not changed since the initial call with
+       @a reuse_inverse false. */
+   void FindInterpolant(Array<Vector*> &x, bool reuse_inverse = false);
 
    /** Set @a diff, comprised of knots in @a kv not contained in this KnotVector.
        @a kv must be of the same order as this KnotVector. The current
@@ -202,6 +207,14 @@ public:
    /** Flag to indicate whether the KnotVector has been coarsened, which means
        it is ready for non-nested refinement. */
    bool coarse;
+
+#ifdef MFEM_USE_LAPACK
+   // Data for reusing banded matrix factorization in FindInterpolant().
+   DenseMatrix fact_AB; /// Banded matrix factorization
+   Array<int> fact_ipiv; /// Row pivot indices
+#else
+   DenseMatrix A_coll_inv; /// Collocation matrix inverse
+#endif
 };
 
 
@@ -285,7 +298,7 @@ public:
        includes the weight. The array of control point coordinates stores each
        point's coordinates contiguously, and points are ordered in a standard
        ijk grid ordering. */
-   NURBSPatch(Array<const KnotVector *> &kv_,  int dim_,
+   NURBSPatch(Array<const KnotVector *> &kv_, int dim_,
               const real_t* control_points);
 
    /// Constructor for a patch of dimension equal to the size of @a kv.
@@ -700,7 +713,8 @@ public:
 
    NURBSExtension(Mesh *mesh_array[], int num_pieces);
 
-   NURBSExtension(const Mesh *patch_topology, const Array<const NURBSPatch*> p);
+   NURBSExtension(const Mesh *patch_topology,
+                  const Array<const NURBSPatch*> &patches_);
 
    /// Copy assignment not supported.
    NURBSExtension& operator=(const NURBSExtension&) = delete;

mesh/submesh/submesh_utils.cpp

@@ -257,15 +257,7 @@ template <typename SubMeshT>
 void AddBoundaryElements(SubMeshT &mesh,
                          const std::unordered_map<int,int> &lface_to_boundary_attribute)
 {
-   mesh.Dimension();
-   const int num_codim_1 = [&mesh]()
-   {
-      auto Dim = mesh.Dimension();
-      if (Dim == 1) { return mesh.GetNV(); }
-      else if (Dim == 2) { return mesh.GetNEdges(); }
-      else if (Dim == 3) { return mesh.GetNFaces(); }
-      else { MFEM_ABORT("Invalid dimension."); return -1; }
-   }();
+   const int num_codim_1 = mesh.GetNumFaces();
 
    if (mesh.Dimension() == 3)
    {

miniapps/nurbs/CMakeLists.txt

@@ -80,6 +80,10 @@ add_mfem_miniapp(nurbs_solenoidal
   LIBRARIES mfem)
   add_dependencies(nurbs_solenoidal copy_miniapps_nurbs_data)
 
+add_mfem_miniapp(nurbs_surface
+  MAIN nurbs_surface.cpp
+  LIBRARIES mfem)
+
 if (MFEM_ENABLE_TESTING)
   add_test(NAME nurbs_ex1_1d_r1_o2_ser
     COMMAND $<TARGET_FILE:nurbs_ex1> -no-vis
@@ -247,6 +251,14 @@ if (MFEM_ENABLE_TESTING)
     COMMAND $<TARGET_FILE:nurbs_solenoidal> -no-vis
         -m ${PROJECT_SOURCE_DIR}/data/cube-nurbs.mesh -r 1 -o 2)
 
+  add_test(NAME nurbs_surface_10_10_10_10_ex1_o3_ser
+    COMMAND $<TARGET_FILE:nurbs_surface> -no-vis
+        -o 3 -nx 10 -ny 10 -fnx 10 -fny 10 -ex 1 -orig)
+
+  add_test(NAME nurbs_surface_10_10_40_40_ex1_o3_ser
+    COMMAND $<TARGET_FILE:nurbs_surface> -no-vis
+        -o 3 -nx 10 -ny 10 -fnx 40 -fny 14 -ex 1)
+
 endif()
 
 if (MFEM_USE_MPI)

miniapps/nurbs/makefile

@@ -21,7 +21,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_MINIAPPS = nurbs_ex1 nurbs_patch_ex1 nurbs_ex3 nurbs_ex5 nurbs_ex24 \
-   nurbs_curveint nurbs_printfunc nurbs_solenoidal nurbs_naca_cmesh
+   nurbs_curveint nurbs_printfunc nurbs_solenoidal nurbs_naca_cmesh nurbs_surface
 PAR_MINIAPPS = nurbs_ex1p nurbs_ex11p
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
@@ -158,6 +158,13 @@ nurbs_naca_cmesh-test-seq: nurbs_naca_cmesh
 nurbs_printfunc-test-seq: nurbs_printfunc
 	@$(call mfem-test,$<,, NURBS miniapp)
 
+SURF_ARGS_1 := -o 3 -nx 10 -ny 10 -fnx 10 -fny 10 -ex 1 -orig
+SURF_ARGS_2 := -o 3 -nx 10 -ny 10 -fnx 40 -fny 40 -ex 1
+
+nurbs_surface-test-seq: nurbs_surface
+	@$(call mfem-test,$<,, NURBS miniapp,$(SURF_ARGS_1))
+	@$(call mfem-test,$<,, NURBS miniapp,$(SURF_ARGS_2))
+
 EX1P_ARGS_1 :=
 EX1P_ARGS_2 := -m ../../data/pipe-nurbs-2d.mesh -o 2 -no-ibp
 EX1P_ARGS_3 := -m ../../data/ball-nurbs.mesh -o 2 --weak-bc -r 0
@@ -192,6 +199,6 @@ clean-build:
 
 clean-exec:
 	@rm -f refined.mesh sin-fit.mesh ex5.mesh exsol.mesh mesh.* sol.* mode_*
-	@rm -f naca-cmesh.mesh sol_?.gf
+	@rm -f naca-cmesh.mesh sol_?.gf *-Surface.mesh
 	@rm -rf Example1* Example3* Example5* Solenoidal_* ParaView
 	@rm -rf CurveInt Naca_cmesh glvis_naca-cmesh.mesh solution.dat

miniapps/nurbs/nurbs_surface.cpp

@@ -0,0 +1,655 @@
+// Copyright (c) 2010-2025, 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.
+//
+//          --------------------------------------------------------
+//          NURBS Surface: Interpolate a 3D Surface in a NURBS Patch
+//          --------------------------------------------------------
+//
+// Compile with: make nurbs_surface
+//
+// Sample runs:  nurbs_surface -o 3 -nx 10 -ny 10 -fnx 10 -fny 10 -ex 1 -orig
+//               nurbs_surface -o 3 -nx 10 -ny 10 -fnx 40 -fny 40 -ex 1
+//               nurbs_surface -o 3 -nx 20 -ny 20 -fnx 10 -fny 10 -ex 1
+//               nurbs_surface -o 3 -nx 20 -ny 20 -fnx 40 -fny 40 -ex 1 -j 0.5
+//               nurbs_surface -o 3 -nx 10 -ny 10 -fnx 10 -fny 10 -ex 2 -orig
+//               nurbs_surface -o 3 -nx 10 -ny 10 -fnx 40 -fny 40 -ex 2
+//               nurbs_surface -o 3 -nx 20 -ny 20 -fnx 10 -fny 10 -ex 2
+//               nurbs_surface -o 3 -nx 10 -ny 10 -fnx 10 -fny 10 -ex 3 -orig
+//               nurbs_surface -o 3 -nx 10 -ny 10 -fnx 40 -fny 40 -ex 3
+//               nurbs_surface -o 3 -nx 20 -ny 20 -fnx 10 -fny 10 -ex 3
+//               nurbs_surface -o 3 -nx 20 -ny 10 -fnx 20 -fny 10 -ex 4 -orig
+//             * nurbs_surface -o 3 -nx 20 -ny 10 -fnx 80 -fny 40 -ex 4
+//             * nurbs_surface -o 3 -nx 40 -ny 20 -fnx 20 -fny 10 -ex 4
+//             * nurbs_surface -o 3 -nx 100 -ny 100 -fnx 100 -fny 100 -ex 5 -orig
+//             * nurbs_surface -o 3 -nx 100 -ny 100 -fnx 400 -fny 400 -ex 5
+//             * nurbs_surface -o 3 -nx 200 -ny 200 -fnx 100 -fny 100 -ex 5
+//
+// Description:  This example demonstrates the use of MFEM to interpolate an
+//               input surface point grid in 3D using a NURBS surface. The NURBS
+//               surface can then be sampled to generate an output mesh of
+//               arbitrary resolution while staying close to the input geometry.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Example data for 3D point grid on surface, given by an analytic function.
+void SurfaceGridExample(int example, int nx, int ny, Array3D<real_t> &vertices,
+                        real_t jitter);
+
+// Write a linear surface mesh with given vertex positions in v.
+void WriteLinearMesh(int nx, int ny, const Array3D<real_t> &v,
+                     const std::string &basename, bool visualization = false,
+                     int x = 0, int y = 0, int w = 500, int h = 500);
+
+// Given an input grid of 3D points on a surface, this class computes a NURBS
+// surface of given order that interpolates the vertices of the input grid.
+class SurfaceInterpolator
+{
+public:
+   /// Constructor for a given 2D point grid size and NURBS order.
+   SurfaceInterpolator(int num_elem_x, int num_elem_y, int order);
+
+   /// Create a surface interpolating the 2D grid of 3D points in @a input3D.
+   void CreateSurface(const Array3D<real_t> &input3D);
+
+   /// Sample the surface with the given grid size, storing points in
+   /// @a output3D.
+   void SampleSurface(int num_elem_x, int num_elem_y, bool compareOriginal,
+                      Array3D<real_t> &output3D);
+
+   /** @brief Write the NURBS surface mesh to file, defined coordinate-wise by
+       the entries of @a cmesh. */
+   void WriteNURBSMesh(const std::string &basename, bool visualization = false,
+                       int x = 0, int y = 0, int w = 500, int h = 500);
+
+protected:
+   /** @brief Compute the NURBS mesh interpolating the given coordinate of the
+       grid of 3D points in @a input3D. */
+   void ComputeNURBS(int coordinate, const Array3D<real_t> &input3D);
+
+private:
+   int nx, ny; // Number of elements in two directions of the surface grid
+   int orderNURBS; // NURBS degree
+   real_t hx, hy, hz; // Grid size in reference space
+
+   Array3D<real_t> initial3D; // Initial grid of points
+
+   static constexpr int dim = 3;
+   Array<int> ncp; // Number of control points in each direction
+   Array<int> nks; // Number of knot-spans in each direction
+
+   std::vector<Vector> ugrid; // Parameter space [0,1]^2 grid point coordinates
+
+   std::vector<KnotVector> kv; // KnotVectors in each direction
+
+   std::unique_ptr<NURBSPatch> patch; // Pointer to the only patch in the mesh
+
+   Mesh mesh; // NURBS mesh representing the surface
+   std::vector<Mesh> cmesh; // NURBS meshes representing point components
+};
+
+
+int main(int argc, char *argv[])
+{
+   // Parse command-line options
+   int nx = 4;
+   int ny = 4;
+   int fnx = 40;
+   int fny = 40;
+   int order = 3;
+   int example = 1;
+   bool visualization = true;
+   bool compareOriginal = false;
+   real_t jitter = 0.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&example, "-ex", "--example",
+                  "Example data");
+   args.AddOption(&nx, "-nx", "--nx",
+                  "Number of elements in x");
+   args.AddOption(&ny, "-ny", "--ny",
+                  "Number of elements in y");
+   args.AddOption(&fnx, "-fnx", "--fnx",
+                  "Number of resampled elements in x");
+   args.AddOption(&fny, "-fny", "--fny",
+                  "Number of resampled elements in y");
+   args.AddOption(&order, "-o", "--order",
+                  "NURBS finite element order (polynomial degree)");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&compareOriginal, "-orig", "--compare-original", "-no-orig",
+                  "--no-compare-original",
+                  "Compare to the original mesh?");
+   args.AddOption(&jitter, "-j", "--jitter",
+                  "Relative jittering in (0,1) to add to the input point "
+                  "coordinates on a uniform nx x ny grid (0 by default)");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   if (compareOriginal && (fnx != nx || fny != ny))
+   {
+      cout << "Comparing to the original mesh requires the same number of "
+           << "samples!\n";
+      return 1;
+   }
+
+   // Dimensions of the 3 surfaces (Input, NURBS, Output)
+   cout << "Input Surface:  " << nx << " x " << ny << " linear elements\n";
+   cout << "NURBS Surface:  " << nx + 1 - order << " x " << ny + 1 - order
+        << " knot elements of order " << order << "\n";
+   cout << "Output Surface: " << fnx << " x " << fny << " linear elements\n";
+
+   // Set the vertex coordinates of the initial linear mesh
+   constexpr int dim = 3;
+   Array3D<real_t> input3D(nx + 1, ny + 1, dim);
+   SurfaceGridExample(example, nx, ny, input3D, jitter);
+
+   // Create a NURBS surface for the given nx, ny and order parameters that
+   // interpolates the input vertex coordinates
+   SurfaceInterpolator surf(nx, ny, order);
+   surf.CreateSurface(input3D);
+
+   // Compute the vertex coordinates of the output linear mesh by sampling the
+   // values from the NURBS surface
+   Array3D<real_t> output3D(fnx + 1, fny + 1, dim);
+   surf.SampleSurface(fnx, fny, compareOriginal, output3D);
+
+   // Save and optionally visualize the 3 surfaces (Input, NURBS, Output)
+   WriteLinearMesh(nx, ny, input3D, "Input-Surface", visualization, 0, 0);
+   surf.WriteNURBSMesh("NURBS-Surface", visualization, 502, 0);
+   WriteLinearMesh(fnx, fny, output3D, "Output-Surface", visualization, 1004, 0);
+
+   return 0;
+}
+
+// f(x,y) = sin(2 * pi * x) * sin(2 * pi * y)
+void Function1(real_t u, real_t v, real_t &x, real_t &y, real_t &z)
+{
+   x = u;
+   y = v;
+   z = sin(2.0 * M_PI * u) * sin(2.0 * M_PI * v);
+}
+
+// Part of the parametric surface of a sphere, using spherical coordinates.
+void Function2(real_t u, real_t v, real_t &x, real_t &y, real_t &z)
+{
+   constexpr real_t r = 1.0;
+   constexpr real_t pi_4 = M_PI * 0.25;
+   constexpr real_t phi0 = -3*pi_4;
+   constexpr real_t phi1 = 3*pi_4;
+   constexpr real_t theta0 = pi_4;
+   constexpr real_t theta1 = 3 * pi_4;
+
+   const real_t phi = (phi0 * (1.0 - v)) + (phi1 * v);
+   const real_t theta = (theta0 * (1.0 - u)) + (theta1 * u);
+   x = r * sin(theta) * cos(phi);
+   y = r * sin(theta) * sin(phi);
+   z = r * cos(theta);
+}
+
+// Helicoid surface
+void Function3(real_t u, real_t v, real_t &x, real_t &y, real_t &z)
+{
+   x = u * cos(2.0 * M_PI * v);
+   y = u * sin(2.0 * M_PI * v);
+   z = v;
+}
+
+// Mobius strip
+void Function4(real_t u, real_t v, real_t &x, real_t &y, real_t &z)
+{
+   constexpr int twists = 1;
+   const real_t a = 1.0 + 0.5 * ((2.0 * v) - 1.0) * cos(2.0 * M_PI * twists * u);
+   x = a * cos(2.0 * M_PI * u);
+   y = a * sin(2.0 * M_PI * u);
+   z = 0.5 * (2.0 * v - 1.0) * sin(2.0 * M_PI * twists * u);
+}
+
+// Breather surface
+void Function5(real_t u, real_t v, real_t &x, real_t &y, real_t &z)
+{
+   const real_t m = 13.2 * ((2.0 * u) - 1.0);
+   const real_t n = 37.4 * ((2.0 * v) - 1.0);
+   constexpr real_t b = 0.4;
+   constexpr real_t r = 1.0 - (b*b);
+   const real_t w = sqrt(r);
+   const real_t denom = b * (pow(w*cosh(b*m),2) + pow(b*sin(w*n),2));
+   x = -m + (2*r*cosh(b*m)*sinh(b*m)) / denom;
+   y = (2*w*cosh(b*m)*(-(w*cos(n)*cos(w*n)) - sin(n)*sin(w*n))) / denom;
+   z = (2*w*cosh(b*m)*(-(w*sin(n)*cos(w*n)) + cos(n)*sin(w*n))) / denom;
+}
+
+void SurfaceFunction(int example, real_t u, real_t v,
+                     real_t &x, real_t &y, real_t &z)
+{
+   switch (example)
+   {
+      case 1:
+         Function1(u, v, x, y, z);
+         break;
+      case 2:
+         Function2(u, v, x, y, z);
+         break;
+      case 3:
+         Function3(u, v, x, y, z);
+         break;
+      case 4:
+         Function4(u, v, x, y, z);
+         break;
+      default:
+         Function5(u, v, x, y, z);
+   };
+}
+
+// Example data for 3D point grid on surface, given by an analytic function.
+void SurfaceExample(int example, const std::vector<Vector> &grid,
+                    Array3D<real_t> &v3D, real_t jitter)
+{
+   int seed = (int)time(0);
+   srand((unsigned)seed);
+
+   real_t h0 = grid[0][1]-grid[0][0], h1 = grid[1][1]-grid[1][0];
+   for (int i = 0; i < grid[0].Size(); i++)
+   {
+      for (int j = 0; j < grid[1].Size(); j++)
+      {
+         if (i != 0 && i != grid[0].Size()-1 && j != 0 && j != grid[1].Size()-1)
+         {
+            SurfaceFunction(example, grid[0][i] + rand_real()*h0*jitter,
+                            grid[1][j] + rand_real()*h1*jitter,
+                            v3D(i, j, 0), v3D(i, j, 1), v3D(i, j, 2));
+         }
+         else
+         {
+            SurfaceFunction(example, grid[0][i], grid[1][j],
+                            v3D(i, j, 0), v3D(i, j, 1), v3D(i, j, 2));
+         }
+      }
+   }
+}
+
+void SurfaceGridExample(int example, int nx, int ny, Array3D<real_t> &vertices,
+                        real_t jitter = 0)
+{
+   // Define a uniform grid of the reference parameter space [0,1]^2
+   std::vector<Vector> uniformGrid(2);
+   for (int i = 0; i < 2; ++i)
+   {
+      const int n = (i == 0) ? nx : ny;
+      const real_t h = 1.0 / n;
+      uniformGrid[i].SetSize(n + 1);
+      for (int j = 0; j <= n; ++j) { uniformGrid[i][j] = j * h; }
+   }
+
+   SurfaceExample(example, uniformGrid, vertices, jitter);
+}
+
+// Write a linear surface mesh with given vertex positions in v.
+void WriteLinearMesh(int nx, int ny, const Array3D<real_t> &v,
+                     const std::string &basename, bool visualization,
+                     int x, int y, int w, int h)
+{
+   const int nv = (nx + 1) * (ny + 1);
+   const int nelem = nx * ny;
+   constexpr int dim = 3; // Spatial dimension
+
+   Mesh lmesh(2, nv, nelem, 0, dim);
+   Vector vertex(dim);
+
+   for (int i = 0; i <= nx; ++i)
+   {
+      for (int j = 0; j <= ny; ++j)
+      {
+         for (int k = 0; k < dim; ++k) { vertex[k] = v(i, j, k); }
+         lmesh.AddVertex(vertex);
+      }
+   }
+
+   Array<int> verts(4);
+
+   auto vID = [&](int i, int j)
+   {
+      return j + (i * (ny + 1));
+   };
+
+   for (int i = 0; i < nx; ++i)
+   {
+      for (int j = 0; j < ny; ++j)
+      {
+         verts[0] = vID(i, j);
+         verts[1] = vID(i+1, j);
+         verts[2] = vID(i+1, j+1);
+         verts[3] = vID(i, j+1);
+
+         Element* el = lmesh.NewElement(Element::QUADRILATERAL);
+         el->SetVertices(verts);
+         lmesh.AddElement(el);
+      }
+   }
+
+   lmesh.FinalizeTopology();
+
+   ofstream mesh_ofs(basename + ".mesh");
+   mesh_ofs.precision(8);
+   lmesh.Print(mesh_ofs);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      constexpr int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "mesh\n" << lmesh
+               << "window_title '" << basename << "'"
+               << "window_geometry "
+               << x << " " << y << " " << w << " " << h << "\n"
+               << "keys PPPPPPPPAattttt******\n"
+               << flush;
+   }
+}
+
+
+// Compute error of interpolation with respect to an input grid of point data.
+void CheckError(const Array3D<real_t> &a, const Array3D<real_t> &b, int c,
+                int nx, int ny)
+{
+   real_t maxErr = 0.0;
+   for (int i = 0; i <= nx; ++i)
+   {
+      for (int j = 0; j <= ny; ++j)
+      {
+         const real_t err_ij = std::abs(a(i, j, c) - b(i, j, 2));
+         maxErr = std::max(maxErr, err_ij);
+      }
+   }
+
+   cout << "Max error: " << maxErr << " for coordinate " << c << endl;
+}
+
+
+// Sample a NURBS mesh to generate a first-order mesh.
+void SampleNURBS(bool uniform, int nx, int ny, const Mesh &mesh,
+                 const Array<int> &nks, const std::vector<Vector> &ugrid,
+                 Array3D<real_t> &vpos)
+{
+   const GridFunction *nodes = mesh.GetNodes();
+
+   const real_t hx = 1.0 / (real_t) nx;
+   const real_t hy = 1.0 / (real_t) ny;
+
+   const real_t hxks = 1.0 / (real_t) nks[0];
+   const real_t hyks = 1.0 / (real_t) nks[1];
+
+   Vector vertex;
+   IntegrationPoint ip;
+
+   ip.z = 1.0;
+   for (int i = 0; i <= nx; ++i)
+   {
+      const real_t xref = uniform ? i * hx : ugrid[0][i];
+      const int nurbsElem0 = std::min((int) (xref / hxks), nks[0] - 1);
+      const real_t ipx = (xref - (nurbsElem0 * hxks)) / hxks;
+      ip.x = ipx;
+
+      for (int j = 0; j <= ny; ++j)
+      {
+         const real_t yref = uniform ? j * hy : ugrid[1][j];
+         const int nurbsElem1 = std::min((int) (yref / hyks), nks[1] - 1);
+         const real_t ipy = (yref - (nurbsElem1 * hyks)) / hyks;
+         ip.y = ipy;
+
+         const int nurbsElem = nurbsElem0 + (nurbsElem1 * nks[0]);
+         nodes->GetVectorValue(nurbsElem, ip, vertex);
+
+         for (int k = 0; k < 3; ++k)
+         {
+            vpos(i, j, k) = vertex[k];
+         }
+      }
+   }
+}
+
+
+SurfaceInterpolator::SurfaceInterpolator(int num_elem_x, int num_elem_y,
+                                         int order) :
+   nx(num_elem_x), ny(num_elem_y), orderNURBS(order),
+   ncp(dim), nks(dim), ugrid(dim - 1)
+{
+   ncp[0] = nx + 1;
+   ncp[1] = ny + 1;
+   ncp[2] = order + 1;
+
+   for (int i = 0; i < dim; ++i)
+   {
+      nks[i] = ncp[i] - order;
+
+      Vector intervals(nks[i]);
+      Array<int> continuity(nks[i] + 1);
+
+      intervals = 1.0 / (real_t) nks[i];
+      continuity = order - 1;
+      continuity[0] = -1;
+      continuity[nks[i]] = -1;
+
+      kv.emplace_back(order, intervals, continuity);
+   }
+
+   patch.reset(new NURBSPatch(&kv[0], &kv[1], &kv[2], dim + 1));
+
+   hx = 1.0 / (real_t) (ncp[0] - 1);
+   hy = 1.0 / (real_t) (ncp[1] - 1);
+   hz = 1.0 / (real_t) (ncp[2] - 1);
+
+   Vector xi_args;
+   Array<int> i_args;
+   for (int i = 0; i < 2; ++i)
+   {
+      kv[i].FindMaxima(i_args, xi_args, ugrid[i]);
+   }
+}
+
+void SurfaceInterpolator::CreateSurface(const Array3D<real_t> &input3D)
+{
+   cmesh.clear();
+   for (int c = 0; c < dim; ++c) // Loop over coordinates
+   {
+      ComputeNURBS(c, input3D);
+      cmesh.emplace_back(mesh);
+   }
+
+   initial3D = input3D;
+}
+
+void SurfaceInterpolator::SampleSurface(int num_elem_x, int num_elem_y,
+                                        bool compareOriginal,
+                                        Array3D<real_t> &output3D)
+{
+   Array3D<real_t> vpos(num_elem_x + 1, num_elem_y + 1, dim);
+   for (int c = 0; c < dim; ++c) // Loop over coordinates
+   {
+      SampleNURBS(true, num_elem_x, num_elem_y, cmesh[c], nks, ugrid, vpos);
+
+      if (compareOriginal)
+      {
+         SampleNURBS(false, num_elem_x, num_elem_y, cmesh[c], nks, ugrid, vpos);
+         CheckError(initial3D, vpos, c, nx, ny);
+      }
+
+      for (int i = 0; i <= num_elem_x; ++i)
+      {
+         for (int j = 0; j <= num_elem_y; ++j)
+         {
+            output3D(i,j,c) = vpos(i,j,2);
+         }
+      }
+   }
+}
+
+void SurfaceInterpolator::ComputeNURBS(int coordinate,
+                                       const Array3D<real_t> &input3D)
+{
+   Array<Vector*> x;
+   for (int i = 0; i < dim; ++i) { x.Append(new Vector(ncp[0])); }
+
+   for (int k = 0; k < ncp[2]; ++k)
+   {
+      const real_t z = k * hz;
+
+      // For each horizontal slice (fixed k), interpolate a 2D surface by
+      // sweeping curve interpolations in each direction. See Algorithm A9.4 of
+      // "The NURBS Book" - 2nd ed - Piegl and Tiller.
+
+      // Resize for sweep in first direction
+      for (int i = 0; i < dim; ++i) { x[i]->SetSize(ncp[0]); }
+
+      // Sweep in the first direction
+      for (int j = 0; j < ncp[1]; ++j)
+      {
+         for (int i = 0; i < ncp[0]; i++)
+         {
+            (*x[0])[i] = ugrid[0][i];
+            (*x[1])[i] = ugrid[1][j];
+
+            const real_t s_ij = input3D(i, j, coordinate);
+            (*x[2])[i] = -1.0 + z + s_ij;
+         }
+
+         const bool reuse_factorization = j > 0;
+         kv[0].FindInterpolant(x, reuse_factorization);
+
+         for (int i = 0; i < ncp[0]; i++)
+         {
+            (*patch)(i,j,k,0) = (*x[0])[i];
+            (*patch)(i,j,k,1) = (*x[1])[i];
+            (*patch)(i,j,k,2) = (*x[2])[i];
+            (*patch)(i,j,k,3) = 1.0; // weight
+         }
+      }
+
+      // Resize for sweep in second direction
+      for (int i = 0; i < dim; ++i) { x[i]->SetSize(ncp[1]); }
+
+      // Do another sweep in the second direction
+      for (int i = 0; i < ncp[0]; i++)
+      {
+         for (int j = 0; j < ncp[1]; ++j)
+         {
+            (*x[0])[j] = (*patch)(i,j,k,0);
+            (*x[1])[j] = (*patch)(i,j,k,1);
+            (*x[2])[j] = (*patch)(i,j,k,2);
+         }
+
+         const bool reuse_factorization = i > 0;
+         kv[1].FindInterpolant(x, reuse_factorization);
+
+         for (int j = 0; j < ncp[1]; ++j)
+         {
+            (*patch)(i,j,k,0) = (*x[0])[j];
+            (*patch)(i,j,k,1) = (*x[1])[j];
+            (*patch)(i,j,k,2) = (*x[2])[j];
+         }
+      }
+   }
+
+   for (auto p : x) { delete p; }
+
+   Array<const NURBSPatch*> patches(1);
+   patches[0] = patch.get();
+   Mesh patch_topology = Mesh::MakeCartesian3D(1, 1, 1, Element::HEXAHEDRON);
+   NURBSExtension nurbsExt(&patch_topology, patches);
+
+   mesh = Mesh(nurbsExt);
+}
+
+void SurfaceInterpolator::WriteNURBSMesh(const std::string &basename,
+                                         bool visualization,
+                                         int x, int y, int w, int h)
+{
+   GridFunction *nodes = cmesh[0].GetNodes();
+   NURBSPatch patch2D(&kv[0], &kv[1], dim);
+   Array<const NURBSPatch*> patches(1);
+   patches[0] = &patch2D;
+   Mesh patch_topology = Mesh::MakeCartesian2D(1, 1, Element::QUADRILATERAL);
+   Array<int> dofs;
+   cmesh[0].NURBSext->GetPatchDofs(0, dofs);
+
+   MFEM_VERIFY(dofs.Size() == (nx + 1) * (ny + 1) * (orderNURBS + 1), "");
+
+   for (int j = 0; j < ncp[1]; ++j)
+   {
+      for (int i = 0; i < ncp[0]; i++)
+      {
+         const int dof = dofs[i + (ncp[0] * (j + (ncp[1] * orderNURBS)))];
+         for (int k = 0; k < 2; ++k) { patch2D(i,j,k) = (*nodes)[dim*dof + k]; }
+         patch2D(i,j,2) = 1.0; // weight
+      }
+   }
+
+   NURBSExtension nurbsExt(&patch_topology, patches);
+   Mesh mesh2D(nurbsExt);
+
+   FiniteElementCollection *fec = nodes->OwnFEC();
+   FiniteElementSpace fespace(&mesh2D, fec, dim, Ordering::byVDIM);
+   GridFunction nodes2D(&fespace);
+
+   const int n = mesh2D.GetNodes()->Size() / (dim - 1);
+   MFEM_VERIFY((dim - 1) * n == mesh2D.GetNodes()->Size(), "");
+   MFEM_VERIFY(dim * n == nodes2D.Size(), "");
+
+   Array<int> dofs2D;
+   mesh2D.NURBSext->GetPatchDofs(0, dofs2D);
+
+   for (int k = 0; k < dim; ++k)
+   {
+      const GridFunction &nodes_k = *cmesh[k].GetNodes();
+
+      for (int j = 0; j < ncp[1]; ++j)
+      {
+         for (int i = 0; i < ncp[0]; i++)
+         {
+            const int dof = dofs[i + (ncp[0] * (j + (ncp[1] * orderNURBS)))];
+            const int dof2D = dofs2D[i + (ncp[0] * j)];
+            nodes2D[(dim*dof2D) + k] = nodes_k[dim*dof + 2];
+         }
+      }
+   }
+
+   // Make mesh2D into a surface mesh with nodes given by nodes2D
+   mesh2D.NewNodes(nodes2D);
+
+   ofstream mesh_ofs(basename + ".mesh");
+   mesh_ofs.precision(8);
+   mesh2D.Print(mesh_ofs);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      constexpr int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "mesh\n" << mesh2D
+               << "window_title '" << basename << "'"
+               << "window_geometry "
+               << x << " " << y << " " << w << " " << h << "\n"
+               << "keys PPPPPPPPAattttt******\n"
+               << flush;
+   }
+}

tests/unit/CMakeLists.txt

@@ -32,6 +32,7 @@ set(UNIT_TESTS_SRCS
   linalg/test_chebyshev.cpp
   linalg/test_complex_dense_matrix.cpp
   linalg/test_complex_operator.cpp
+  linalg/test_complex_vector.cpp
   linalg/test_constrainedsolver.cpp
   linalg/test_direct_solvers.cpp
   linalg/test_hypre_ilu.cpp

tests/unit/linalg/test_complex_vector.cpp

@@ -0,0 +1,145 @@
+// Copyright (c) 2010-2025, 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 "mfem.hpp"
+#include "unit_tests.hpp"
+#include <numeric>
+
+using namespace mfem;
+
+TEST_CASE("Complex Vector init-list and C-style array constructors",
+          "[ComplexVector]")
+{
+   std::complex<real_t> ContigData[6] = {std::complex<real_t>(6.0,1.0),
+                                         std::complex<real_t>(5.0,2.0),
+                                         std::complex<real_t>(4.0,3.0),
+                                         std::complex<real_t>(3.0,4.0),
+                                         std::complex<real_t>(2.0,5.0),
+                                         std::complex<real_t>(1.0,6.0)
+                                        };
+   // Point and size constructor
+   ComplexVector a(ContigData, 6);
+   // Braced-list constructor
+   ComplexVector b({std::complex<real_t>(6.0,1.0),
+                    std::complex<real_t>(5.0,2.0),
+                    std::complex<real_t>(4.0,3.0),
+                    std::complex<real_t>(3.0,4.0),
+                    std::complex<real_t>(2.0,5.0),
+                    std::complex<real_t>(1.0,6.0)});
+   // Statically sized C-style array constructor
+   ComplexVector c(ContigData);
+
+   for (int i = 0; i < a.Size(); i++)
+   {
+      REQUIRE(a[i] == b[i]);
+      REQUIRE(a[i] == c[i]);
+   }
+}
+
+TEST_CASE("Complex Vector Move Constructor", "[ComplexVector]")
+{
+   constexpr int N = 6;
+   std::complex<real_t> ContigData[6] = {std::complex<real_t>(6.0,1.0),
+                                         std::complex<real_t>(5.0,2.0),
+                                         std::complex<real_t>(4.0,3.0),
+                                         std::complex<real_t>(3.0,4.0),
+                                         std::complex<real_t>(2.0,5.0),
+                                         std::complex<real_t>(1.0,6.0)
+                                        };
+   ComplexVector a(ContigData, N);
+   ComplexVector b(N);
+   for (int i = 0; i < N; i++)
+   {
+      b(i) = std::complex<real_t>(N - i, i + 1);
+   }
+
+   std::complex<real_t>* a_data = a.GetData();
+   std::complex<real_t>* b_data = b.GetData();
+
+   ComplexVector move_non_owning(std::move(a));
+   ComplexVector move_owning(std::move(b));
+
+   REQUIRE(a.Size() == 0);
+   REQUIRE(a.GetData() == nullptr);
+   REQUIRE(b.Size() == 0);
+   REQUIRE(b.GetData() == nullptr);
+
+   // Should both be no-ops
+   a.Destroy();
+   b.Destroy();
+
+   REQUIRE(move_non_owning.OwnsData() == false);
+   REQUIRE(move_owning.OwnsData() == true);
+
+   REQUIRE(move_non_owning.Size() == N);
+   REQUIRE(move_owning.Size() == N);
+
+   // Make sure that the pointers were reused
+   REQUIRE(move_non_owning.GetData() == a_data);
+   REQUIRE(move_owning.GetData() == b_data);
+
+   for (int i = 0; i < N; i++)
+   {
+      REQUIRE(move_non_owning(i) == std::complex<real_t>(N - i, i + 1));
+      REQUIRE(move_owning(i) == std::complex<real_t>(N - i, i + 1));
+   }
+}
+
+TEST_CASE("Complex Vector Move Assignment", "[ComplexVector]")
+{
+   constexpr int N = 6;
+   std::complex<real_t> ContigData[6] = {std::complex<real_t>(6.0,1.0),
+                                         std::complex<real_t>(5.0,2.0),
+                                         std::complex<real_t>(4.0,3.0),
+                                         std::complex<real_t>(3.0,4.0),
+                                         std::complex<real_t>(2.0,5.0),
+                                         std::complex<real_t>(1.0,6.0)
+                                        };
+   ComplexVector a(ContigData, N);
+   ComplexVector b(N);
+   for (int i = 0; i < N; i++)
+   {
+      b(i) = std::complex<real_t>(N - i, i + 1);
+   }
+
+   std::complex<real_t>* a_data = a.GetData();
+   std::complex<real_t>* b_data = b.GetData();
+
+   ComplexVector move_non_owning;
+   move_non_owning = std::move(a);
+   ComplexVector move_owning;
+   move_owning = std::move(b);
+
+   REQUIRE(a.Size() == 0);
+   REQUIRE(a.GetData() == nullptr);
+   REQUIRE(b.Size() == 0);
+   REQUIRE(b.GetData() == nullptr);
+
+   // Should both be no-ops
+   a.Destroy();
+   b.Destroy();
+
+   REQUIRE(move_non_owning.OwnsData() == false);
+   REQUIRE(move_owning.OwnsData() == true);
+
+   REQUIRE(move_non_owning.Size() == N);
+   REQUIRE(move_owning.Size() == N);
+
+   // Make sure that the pointers were reused
+   REQUIRE(move_non_owning.GetData() == a_data);
+   REQUIRE(move_owning.GetData() == b_data);
+
+   for (int i = 0; i < N; i++)
+   {
+      REQUIRE(move_non_owning(i) == std::complex<real_t>(N - i, i + 1));
+      REQUIRE(move_owning(i) == std::complex<real_t>(N - i, i + 1));
+   }
+}