diff --git a/CHANGELOG b/CHANGELOG
index 81750882..f00453a2 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -23,6 +23,13 @@ Unlike previous GLVis releases, this version requires a C++17 compiler.
 - Added non-persistent mode of the server, when the server terminates after all
   visualization windows are closed.
 
+- Added visualization of complex function data originating from
+  '(Par)ComplexGridFunction' class in MFEM. The complex values are represented
+  by their magnitude, phase, real part or imaginary part, which can be switched
+  by 'Q' key. Additionally, the complex function can be animated by an added
+  harmonic phase. To increase or decrease the phase frequency, use ';' or ':',
+  respectively.
+
 - Implemented DOF numbering in 2D scalar mode, pressing the `n` or `N` key cycles
   through: None → Elements → Edges → Vertices → DOFs. Parallel numbering is now by
   default 'local' to each rank; 'global' vs. 'local' numbering can be toggled with
diff --git a/README.md b/README.md
index afb0016c..ec95b1ea 100644
--- a/README.md
+++ b/README.md
@@ -154,6 +154,11 @@ Key commands
   - piece-wise constant refined (LOR)
   - L2 element dof collocation (interpolation)
   - L2 element projection (L2 projection)
+- <kbd>Q</kbd> – Cycle between representations of the visualized *complex data*. The options are:
+  - magnitude
+  - phase
+  - real part
+  - imaginary part
 - <kbd>q</kbd> – Exit
 
 ## Advanced
@@ -196,6 +201,7 @@ Key commands
   - `on` (default): recompute both the bounding box and the value range
   - `value`: recompute only the value range
   - `mesh`: recompute only the bounding box
+- <kbd>Alt</kbd> + <kbd>.</kbd> – Start/stop harmonic complex phase animation (speed can be controlled with <kbd>Alt</kbd> + <kbd>0</kbd> / <kbd>Alt</kbd> + <kbd>Enter</kbd>)
 
 ## 2D scalar data
 
diff --git a/lib/aux_vis.cpp b/lib/aux_vis.cpp
index 204ad80d..60f77a0b 100644
--- a/lib/aux_vis.cpp
+++ b/lib/aux_vis.cpp
@@ -204,8 +204,8 @@ GLWindow* InitVisualization(const char name[], int x, int y, int w, int h,
    wnd->setOnKeyDown (SDLK_s, KeyS);
    wnd->setOnKeyDown ('S', KeyS);
 
-   wnd->setOnKeyDown (SDLK_q, KeyQPressed);
-   // wnd->setOnKeyDown (SDLK_Q, KeyQPressed);
+   wnd->setOnKeyDown (SDLK_q, KeyqPressed);
+   wnd->setOnKeyDown ('Q', KeyQPressed);
 
    wnd->setOnKeyDown (SDLK_LEFT, KeyLeftPressed);
    wnd->setOnKeyDown (SDLK_RIGHT, KeyRightPressed);
@@ -609,12 +609,28 @@ void RemoveIdleFunc(void (*Func)(void))
 
 
 thread_local double xang = 0., yang = 0.;
+const double xang_step = 0.2; // angle in degrees
 thread_local gl3::GlMatrix srot;
 thread_local double sph_t, sph_u;
 thread_local GLint oldx, oldy, startx, starty;
 
 thread_local int constrained_spinning = 0;
 
+thread_local double phase_rate = 0.;
+thread_local int phase_anim = 0;
+const double phase_step = 0.001;
+
+void CheckMainIdleFunc()
+{
+   if (locscene->spinning || phase_anim)
+   {
+      AddIdleFunc(MainLoop);
+   }
+   if (!locscene->spinning && !phase_anim)
+   {
+      RemoveIdleFunc(MainLoop);
+   }
+}
 
 void MainLoop()
 {
@@ -624,13 +640,19 @@ void MainLoop()
       if (!constrained_spinning)
       {
          locscene->Rotate(xang, yang);
-         SendExposeEvent();
       }
       else
       {
          locscene->PreRotate(xang, 0.0, 0.0, 1.0);
-         SendExposeEvent();
       }
+   }
+   if (phase_anim)
+   {
+      locscene->UpdateComplexPhase(phase_rate);
+   }
+   if (locscene->spinning || phase_anim)
+   {
+      SendExposeEvent();
       std::this_thread::sleep_for(std::chrono::milliseconds{10}); // sleep for 0.01 seconds
    }
    if (locscene->movie)
@@ -674,7 +696,7 @@ inline void ComputeSphereAngles(int &newx, int &newy,
 void LeftButtonDown(GLWindow::MouseEventInfo *event)
 {
    locscene -> spinning = 0;
-   RemoveIdleFunc(MainLoop);
+   CheckMainIdleFunc();
 
    oldx = event->mouse_x;
    oldy = event->mouse_y;
@@ -1334,12 +1356,17 @@ void KeyCtrlP()
 #endif
 }
 
-void KeyQPressed()
+void KeyqPressed()
 {
    wnd->signalQuit();
    visualize = 0;
 }
 
+void KeyQPressed()
+{
+   Window::SwitchSolution();
+}
+
 void ToggleThreads()
 {
    static const char *state[] = { "running", "stopped" };
@@ -1350,7 +1377,7 @@ void ToggleThreads()
    }
 }
 
-void ThreadsPauseFunc(GLenum state)
+void ThreadsPauseFunc(SDL_Keymod state)
 {
    if (state & KMOD_CTRL)
    {
@@ -1386,57 +1413,96 @@ void CheckSpin()
    {
       xang = 0.;
    }
-   if (xang != 0. || yang != 0.)
-   {
-      locscene->spinning = 1;
-      AddIdleFunc(MainLoop);
-   }
-   else
-   {
-      locscene->spinning = 0;
-      RemoveIdleFunc(MainLoop);
-   }
+   locscene->spinning = (xang != 0. || yang != 0.);
+   CheckMainIdleFunc();
    cout << "Spin angle: " << xang << " degrees / frame" << endl;
 }
 
-const double xang_step = 0.2; // angle in degrees
+void CheckPhaseAnim()
+{
+   if (fabs(phase_rate) < phase_step / 2.)
+   {
+      phase_rate = 0.;
+   }
+   phase_anim = (phase_rate != 0.);
+   CheckMainIdleFunc();
+   cout << "Phase rate: " << phase_rate << " period / frame" << endl;
+}
 
-void Key0Pressed()
+void Key0Pressed(SDL_Keymod mod)
 {
-   if (!locscene -> spinning)
+   if (mod & KMOD_ALT)
+   {
+      if (!phase_anim)
+      {
+         phase_rate = 0.;
+      }
+      phase_rate -= phase_step;
+      CheckPhaseAnim();
+   }
+   else
    {
-      xang = 0;
+      if (!locscene -> spinning)
+      {
+         xang = 0;
+      }
+      xang -= xang_step;
+      CheckSpin();
    }
-   xang -= xang_step;
-   CheckSpin();
 }
 
-void KeyDeletePressed()
+void KeyDeletePressed(SDL_Keymod mod)
 {
-   if (locscene -> spinning)
+   if (mod & KMOD_ALT)
    {
-      xang = yang = 0.;
-      locscene -> spinning = 0;
-      RemoveIdleFunc(MainLoop);
-      constrained_spinning = 1;
+      if (phase_anim)
+      {
+         phase_rate = 0.;
+         phase_anim = 0;
+      }
+      else
+      {
+         phase_rate = phase_step;
+         phase_anim = 1;
+      }
    }
    else
    {
-      xang = xang_step;
-      locscene -> spinning = 1;
-      AddIdleFunc(MainLoop);
+      if (locscene -> spinning)
+      {
+         xang = yang = 0.;
+         locscene -> spinning = 0;
+      }
+      else
+      {
+         xang = xang_step;
+         locscene -> spinning = 1;
+      }
       constrained_spinning = 1;
    }
+   CheckMainIdleFunc();
 }
 
-void KeyEnterPressed()
+void KeyEnterPressed(SDL_Keymod mod)
 {
-   if (!locscene -> spinning)
+   if (mod & KMOD_ALT)
    {
-      xang = 0;
+      if (!phase_anim)
+      {
+         phase_rate = 0.;
+      }
+      phase_rate += phase_step;
+      CheckPhaseAnim();
+   }
+   else
+   {
+      if (!locscene -> spinning)
+      {
+         xang = 0;
+      }
+      xang += xang_step;
+      CheckSpin();
    }
-   xang += xang_step;
-   CheckSpin();
 }
 
 void Key7Pressed()
@@ -1515,7 +1581,7 @@ void ShiftView(double dx, double dy)
    locscene->ViewCenterY += dy/scale;
 }
 
-void KeyLeftPressed(GLenum state)
+void KeyLeftPressed(SDL_Keymod state)
 {
    if (state & KMOD_CTRL)
    {
@@ -1528,7 +1594,7 @@ void KeyLeftPressed(GLenum state)
    SendExposeEvent();
 }
 
-void KeyRightPressed(GLenum state)
+void KeyRightPressed(SDL_Keymod state)
 {
    if (state & KMOD_CTRL)
    {
@@ -1541,7 +1607,7 @@ void KeyRightPressed(GLenum state)
    SendExposeEvent();
 }
 
-void KeyUpPressed(GLenum state)
+void KeyUpPressed(SDL_Keymod state)
 {
    if (state & KMOD_CTRL)
    {
@@ -1554,7 +1620,7 @@ void KeyUpPressed(GLenum state)
    SendExposeEvent();
 }
 
-void KeyDownPressed(GLenum state)
+void KeyDownPressed(SDL_Keymod state)
 {
    if (state & KMOD_CTRL)
    {
diff --git a/lib/aux_vis.hpp b/lib/aux_vis.hpp
index f64e1734..ce7ee31e 100644
--- a/lib/aux_vis.hpp
+++ b/lib/aux_vis.hpp
@@ -52,6 +52,7 @@ void SetLegacyGLOnly(bool status);
 
 void AddIdleFunc(void (*Func)(void));
 void RemoveIdleFunc(void (*Func)(void));
+void CheckMainIdleFunc();
 
 void LeftButtonDown  (GLWindow::MouseEventInfo *event);
 void LeftButtonLoc   (GLWindow::MouseEventInfo *event);
@@ -67,9 +68,10 @@ void TouchPinch(SDL_MultiGestureEvent & e);
 
 void KeyCtrlP();
 void KeyS();
+void KeyqPressed();
 void KeyQPressed();
 void ToggleThreads();
-void ThreadsPauseFunc(GLenum);
+void ThreadsPauseFunc(SDL_Keymod);
 void ThreadsStop();
 void ThreadsRun();
 
@@ -83,14 +85,14 @@ void Key7Pressed();
 void Key8Pressed();
 void Key9Pressed();
 
-void Key0Pressed();
-void KeyDeletePressed();
-void KeyEnterPressed();
+void Key0Pressed(SDL_Keymod);
+void KeyDeletePressed(SDL_Keymod);
+void KeyEnterPressed(SDL_Keymod);
 
-void KeyLeftPressed(GLenum);
-void KeyRightPressed(GLenum);
-void KeyUpPressed(GLenum);
-void KeyDownPressed(GLenum);
+void KeyLeftPressed(SDL_Keymod);
+void KeyRightPressed(SDL_Keymod);
+void KeyUpPressed(SDL_Keymod);
+void KeyDownPressed(SDL_Keymod);
 void KeyJPressed();
 void KeyMinusPressed();
 void KeyPlusPressed();
diff --git a/lib/data_state.cpp b/lib/data_state.cpp
index 6fbf1af8..36d2fedd 100644
--- a/lib/data_state.cpp
+++ b/lib/data_state.cpp
@@ -41,12 +41,14 @@ DataState &DataState::operator=(DataState &&ss)
    internal = std::move(ss.internal);
 
    type = ss.type;
+   cmplx_sol = ss.cmplx_sol;
    quad_sol = ss.quad_sol;
    keys = std::move(ss.keys);
 
    fix_elem_orient = ss.fix_elem_orient;
    save_coloring = ss.save_coloring;
    keep_attr = ss.keep_attr;
+   cmplx_phase = ss.cmplx_phase;
 
    return *this;
 }
@@ -62,6 +64,7 @@ void DataState::SetMesh(std::unique_ptr<Mesh> &&pmesh)
    internal.mesh = std::move(pmesh);
    internal.mesh_quad.reset();
    if (grid_f && grid_f->FESpace()->GetMesh() != mesh.get()) { SetGridFunction(NULL); }
+   if (cgrid_f && cgrid_f->FESpace()->GetMesh() != mesh.get()) { SetCmplxGridFunction(NULL); }
    if (quad_f && quad_f->GetSpace()->GetMesh() != mesh.get()) { SetQuadFunction(NULL); }
 }
 
@@ -146,6 +149,8 @@ void DataState::SetGridFunction(std::unique_ptr<GridFunction> &&pgf,
                                 int component)
 {
    internal.grid_f = std::move(pgf);
+   internal.cgrid_f.reset();
+   cmplx_sol = ComplexSolution::NONE;
    internal.quad_f.reset();
    quad_sol = QuadSolution::NONE;
    SetGridFunctionSolution(component);
@@ -156,7 +161,75 @@ void DataState::SetGridFunction(std::vector<GridFunction*> &gf_array,
 {
    SetGridFunction(new GridFunction(mesh.get(), gf_array.data(), num_pieces),
                    component);
-   if (!keep_attr) { ComputeDofsOffsets(gf_array); }
+   if (!keep_attr)
+   {
+      std::vector<const FiniteElementSpace*> fespaces(gf_array.size());
+      for (size_t i = 0; i < gf_array.size(); i++)
+      {
+         fespaces[i] = gf_array[i]->FESpace();
+      }
+      ComputeDofsOffsets(fespaces);
+   }
+}
+
+void DataState::SetCmplxGridFunction(ComplexGridFunction *gf,
+                                     int component)
+{
+   if (cgrid_f.get() != gf)
+   {
+      internal.grid_f.reset();
+      cmplx_sol = ComplexSolution::NONE;
+   }
+   internal.cgrid_f.reset(gf);
+   internal.quad_f.reset();
+   quad_sol = QuadSolution::NONE;
+   SetComplexFunctionSolution(component);
+}
+
+void DataState::SetCmplxGridFunction(
+   std::unique_ptr<ComplexGridFunction> &&pgf, int component)
+{
+   if (cgrid_f.get() != pgf.get())
+   {
+      internal.grid_f.reset();
+      cmplx_sol = ComplexSolution::NONE;
+   }
+   internal.cgrid_f = std::move(pgf);
+   internal.quad_f.reset();
+   quad_sol = QuadSolution::NONE;
+   SetComplexFunctionSolution(component);
+}
+
+void DataState::SetCmplxGridFunction(
+   const std::vector<ComplexGridFunction *> &cgf_array, int component)
+{
+   const int np = cgf_array.size();
+   std::vector<GridFunction *> r_array(np), i_array(np);
+   for (int p = 0; p < np; p++)
+   {
+      r_array[p] = &(cgf_array[p]->real());
+      i_array[p] = &(cgf_array[p]->imag());
+   }
+   GridFunction *rgf = new GridFunction(mesh.get(), r_array.data(), np);
+   GridFunction *igf = new GridFunction(mesh.get(), i_array.data(), np);
+   ComplexGridFunction *cgf = new ComplexGridFunction(rgf->FESpace());
+   // transfer ownership of the FES
+   cgf->MakeOwner(rgf->OwnFEC());
+   rgf->MakeOwner(NULL);
+   cgf->real() = *rgf;
+   cgf->imag() = *igf;
+   delete rgf;
+   delete igf;
+   SetCmplxGridFunction(cgf, component);
+   if (!keep_attr)
+   {
+      std::vector<const FiniteElementSpace*> fespaces(cgf_array.size());
+      for (size_t i = 0; i < cgf_array.size(); i++)
+      {
+         fespaces[i] = cgf_array[i]->FESpace();
+      }
+      ComputeDofsOffsets(fespaces);
+   }
 }
 
 void DataState::SetQuadFunction(QuadratureFunction *qf, int component)
@@ -167,6 +240,7 @@ void DataState::SetQuadFunction(QuadratureFunction *qf, int component)
       quad_sol = QuadSolution::NONE;
    }
    internal.quad_f.reset(qf);
+   internal.cgrid_f.reset();
    SetQuadFunctionSolution(component);
 }
 
@@ -179,6 +253,7 @@ void DataState::SetQuadFunction(std::unique_ptr<QuadratureFunction> &&pqf,
       quad_sol = QuadSolution::NONE;
    }
    internal.quad_f = std::move(pqf);
+   internal.cgrid_f.reset();
    SetQuadFunctionSolution(component);
 }
 
@@ -275,7 +350,29 @@ void DataState::Extrude1DMeshAndSolution()
 
    Mesh *mesh2d = Extrude1D(mesh.get(), 1, 0.1*(xmax - xmin));
 
-   if (grid_f)
+   if (cgrid_f)
+   {
+      if (cgrid_f->VectorDim() > 1)
+      {
+         ProjectVectorFEGridFunction();
+      }
+
+      FiniteElementCollection *fec2d = new L2_FECollection(
+         cgrid_f->FESpace()->FEColl()->GetOrder(), 2);
+      FiniteElementSpace *fes2d = new FiniteElementSpace(mesh2d, fec2d,
+                                                         cgrid_f->FESpace()->GetVDim());
+      ComplexGridFunction *cgrid_f_2d = new ComplexGridFunction(fes2d);
+      cgrid_f_2d->MakeOwner(fec2d);
+      VectorGridFunctionCoefficient vcsol_r(&cgrid_f->real());
+      VectorGridFunctionCoefficient vcsol_i(&cgrid_f->imag());
+      ::VectorExtrudeCoefficient vc2d_r(mesh.get(), vcsol_r, 1);
+      ::VectorExtrudeCoefficient vc2d_i(mesh.get(), vcsol_i, 1);
+      cgrid_f_2d->real().ProjectCoefficient(vc2d_r);
+      cgrid_f_2d->imag().ProjectCoefficient(vc2d_i);
+
+      internal.cgrid_f.reset(cgrid_f_2d);
+   }
+   else if (grid_f)
    {
       if (grid_f->VectorDim() > 1)
       {
@@ -319,7 +416,9 @@ void DataState::Extrude1DMeshAndSolution()
 
 void DataState::Extrude2D3VMeshAndSolution()
 {
-   if (mesh->SpaceDimension() == 3 || !grid_f || grid_f->VectorDim() < 3) { return; }
+   if (mesh->SpaceDimension() == 3) { return; }
+   if (cgrid_f && cgrid_f->VectorDim() < 3) { return; }
+   if (!grid_f || grid_f->VectorDim() < 3) { return; }
 
    // not all vector elements can be embedded in 3D, so conversion to L2 elements
    // is performed already here
@@ -339,14 +438,25 @@ void DataState::Extrude2D3VMeshAndSolution()
       mesh3d->SetCurvature(1, false, 3);
    }
 
-   FiniteElementSpace *fes2d = grid_f->FESpace();
+   FiniteElementSpace *fes2d = (cgrid_f)?(cgrid_f->FESpace()):(grid_f->FESpace());
    FiniteElementSpace *fes3d = new FiniteElementSpace(*fes2d, mesh3d);
-   GridFunction *gf3d = new GridFunction(fes3d);
-   *gf3d = *grid_f;
-   grid_f->MakeOwner(NULL);
-   gf3d->MakeOwner(const_cast<FiniteElementCollection*>(fes2d->FEColl()));
+   if (cgrid_f)
+   {
+      ComplexGridFunction *cgf3d = new ComplexGridFunction(fes3d);
+      (Vector&)*cgf3d = (const Vector&)*cgrid_f;
+      cgrid_f->MakeOwner(NULL);
+      cgf3d->MakeOwner(const_cast<FiniteElementCollection*>(fes2d->FEColl()));
+      SetCmplxGridFunction(cgf3d);
+   }
+   else
+   {
+      GridFunction *gf3d = new GridFunction(fes3d);
+      *gf3d = *grid_f;
+      grid_f->MakeOwner(NULL);
+      gf3d->MakeOwner(const_cast<FiniteElementCollection*>(fes2d->FEColl()));
+      SetGridFunction(gf3d);
+   }
 
-   SetGridFunction(gf3d);
    delete fes2d;
    SetMesh(mesh3d);
 }
@@ -434,6 +544,48 @@ void DataState::SetGridFunctionSolution(int gf_component)
    type = (grid_f->VectorDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR;
 }
 
+void DataState::SetComplexFunctionSolution(int gf_component)
+{
+   if (!cgrid_f)
+   {
+      type = (mesh)?(FieldType::MESH):(FieldType::UNKNOWN);
+      return;
+   }
+
+   if (gf_component != -1)
+   {
+      if (gf_component < 0 || gf_component >= cgrid_f->FESpace()->GetVDim())
+      {
+         cerr << "Invalid component " << gf_component << '.' << endl;
+         exit(1);
+      }
+      FiniteElementSpace *ofes = cgrid_f->FESpace();
+      FiniteElementCollection *fec =
+         FiniteElementCollection::New(ofes->FEColl()->Name());
+      FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec);
+      ComplexGridFunction *new_gf = new ComplexGridFunction(fes);
+      new_gf->MakeOwner(fec);
+      for (int i = 0; i < new_gf->real().Size(); i++)
+      {
+         (new_gf->real())(i) = (cgrid_f->real())(ofes->DofToVDof(i, gf_component));
+         (new_gf->imag())(i) = (cgrid_f->imag())(ofes->DofToVDof(i, gf_component));
+      }
+      SetCmplxGridFunction(new_gf);
+      return;
+   }
+
+   if (cgrid_f->VectorDim() == 1)
+   {
+      type = FieldType::SCALAR;
+   }
+   else
+   {
+      type = FieldType::VECTOR;
+   }
+
+   SetComplexSolution();
+}
+
 void DataState::SetQuadFunctionSolution(int qf_component)
 {
    if (!quad_f)
@@ -474,6 +626,200 @@ void DataState::SetQuadFunctionSolution(int qf_component)
    SetQuadSolution();
 }
 
+void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print)
+{
+   double (*cmplx_2_scalar)(double, double);
+   const char *str_cmplx_2_scalar;
+   bool linear;
+
+   switch (cmplx_type)
+   {
+      case ComplexSolution::Magnitude:
+         cmplx_2_scalar = [](double r, double i) { return hypot(r,i); };
+         str_cmplx_2_scalar = "magnitude";
+         linear = false;
+         break;
+      case ComplexSolution::Phase:
+         cmplx_2_scalar = [](double r, double i) { return atan2(r,i); };
+         str_cmplx_2_scalar = "phase";
+         linear = false;
+         break;
+      case ComplexSolution::Real:
+         cmplx_2_scalar = [](double r, double i) { return r; };
+         str_cmplx_2_scalar = "real part";
+         linear = true;
+         break;
+      case ComplexSolution::Imag:
+         cmplx_2_scalar = [](double r, double i) { return i; };
+         str_cmplx_2_scalar = "imaginary part";
+         linear = true;
+         break;
+      default:
+         cout << "Unknown complex data representation" << endl;
+         return;
+   }
+
+   if (print)
+   {
+      cout << "Representing complex function by: " << str_cmplx_2_scalar;
+      if (cmplx_type != ComplexSolution::Magnitude) { cout << " (+ animated harmonic phase)"; }
+      cout << endl;
+   }
+
+   FiniteElementSpace *fes = cgrid_f->FESpace();
+   const int dim = mesh->Dimension();
+   const int vdim = fes->GetVectorDim();
+
+   const bool scal_fes = (fes->FEColl()->GetMapType(dim) == FiniteElement::VALUE
+                          && fes->FEColl()->GetRangeType(dim) == FiniteElement::SCALAR);
+
+   const double cos_ph = cos(2. * M_PI * cmplx_phase);
+   const double sin_ph = sin(2. * M_PI * cmplx_phase);
+   auto rot_ph = [cos_ph, sin_ph](double &r, double &i)
+   {
+      double rph = +r * cos_ph - i * sin_ph;
+      double iph = +r * sin_ph + i * cos_ph;
+      r = rph, i = iph;
+   };
+
+   GridFunction *gf{};
+   if (scal_fes || linear)
+   {
+      gf = new GridFunction(fes);
+   }
+   else
+   {
+      FiniteElementCollection *l2_fec = new L2_FECollection(fes->FEColl()->GetOrder(),
+                                                            dim);
+      FiniteElementSpace *l2_fes = new FiniteElementSpace(mesh.get(), l2_fec, vdim);
+      gf = new GridFunction(l2_fes);
+      gf->MakeOwner(l2_fec);
+   }
+
+   if (scal_fes)
+   {
+      for (int i = 0; i < fes->GetVSize(); i++)
+      {
+         double rval = cgrid_f->real()(i);
+         double ival = cgrid_f->imag()(i);
+         rot_ph(rval, ival);
+         (*gf)(i) = cmplx_2_scalar(rval, ival);
+      }
+   }
+   else
+   {
+      Array<int> vdofs;
+      Vector r_data, i_data, z_data;
+      DenseMatrix r_vals, i_vals, z_vals;
+
+      if (linear)
+      {
+         // pointwise projection at DOFs
+         DenseMatrix vshape;
+         Vector shape;
+
+         auto kernel = [&](double rdof, double idof, double w)
+         {
+            double rval = rdof * w;
+            double ival = idof * w;
+            rot_ph(rval, ival);
+            const double zval = cmplx_2_scalar(rval, ival);
+            return (w != 0.)?(zval / w):(0.);
+         };
+
+         for (int z = 0; z < mesh->GetNE(); z++)
+         {
+            fes->GetElementVDofs(z, vdofs);
+            cgrid_f->real().GetSubVector(vdofs, r_data);
+            cgrid_f->imag().GetSubVector(vdofs, i_data);
+            z_data.SetSize(vdofs.Size());
+
+            ElementTransformation *Tr = fes->GetElementTransformation(z);
+            const FiniteElement *fe = fes->GetFE(z);
+            const IntegrationRule &ir = fe->GetNodes();
+            const int nnp = ir.GetNPoints();
+
+            if (fe->GetRangeType() == FiniteElement::SCALAR)
+            {
+               shape.SetSize(nnp);
+               r_vals.Reset(r_data.GetData(), nnp, vdim);
+               i_vals.Reset(i_data.GetData(), nnp, vdim);
+               z_vals.Reset(z_data.GetData(), nnp, vdim);
+
+               for (int n = 0; n < nnp; n++)
+               {
+                  const IntegrationPoint &ip = ir.IntPoint(n);
+                  Tr->SetIntPoint(&ip);
+                  fe->CalcPhysShape(*Tr, shape);
+                  const double w = shape(n);
+                  for (int d = 0; d < vdim; d++)
+                  {
+                     z_vals(n,d) = kernel(r_vals(n,d), i_vals(n,d), w);
+                  }
+               }
+            }
+            else
+            {
+               vshape.SetSize(nnp, vdim);
+
+               for (int n = 0; n < nnp; n++)
+               {
+                  const IntegrationPoint &ip = ir.IntPoint(n);
+                  Tr->SetIntPoint(&ip);
+                  fe->CalcPhysVShape(*Tr, vshape);
+                  Vector vec(vdim);
+                  vshape.GetRow(n, vec);
+                  const double w = vec.Norml2();
+                  z_data(n) = kernel(r_data(n), i_data(n), w);
+               }
+            }
+
+            gf->SetSubVector(vdofs, z_data);
+         }
+      }
+      else
+      {
+         // pointwise projection to L2
+
+         const FiniteElementSpace *l2_fes = gf->FESpace();
+         Vector r_vec, i_vec;
+
+         for (int z = 0; z < mesh->GetNE(); z++)
+         {
+            l2_fes->GetElementVDofs(z, vdofs);
+
+            ElementTransformation *Tr = fes->GetElementTransformation(z);
+            const FiniteElement *fe = l2_fes->GetFE(z);
+            const IntegrationRule &ir = fe->GetNodes();
+            const int nnp = ir.GetNPoints();
+            cgrid_f->real().GetVectorValues(*Tr, ir, r_vals);
+            cgrid_f->imag().GetVectorValues(*Tr, ir, i_vals);
+            z_data.SetSize(nnp * vdim);
+            z_vals.Reset(z_data.GetData(), nnp, vdim);
+
+            for (int n = 0; n < nnp; n++)
+            {
+               r_vals.GetColumnReference(n, r_vec);
+               i_vals.GetColumnReference(n, i_vec);
+
+               for (int d = 0; d < vdim; d++)
+               {
+                  double rval = r_vec(d);
+                  double ival = i_vec(d);
+                  rot_ph(rval, ival);
+                  z_vals(n,d) = cmplx_2_scalar(rval, ival);
+               }
+            }
+
+            gf->SetSubVector(vdofs, z_data);
+         }
+      }
+   }
+   internal.grid_f.reset(gf);
+
+   cmplx_sol = cmplx_type;
+}
+
 void DataState::SetQuadSolution(QuadSolution quad_type)
 {
    // assume identical order
@@ -626,13 +972,313 @@ DataState::ProjectVectorFEGridFunction(std::unique_ptr<GridFunction> gf)
    return gf;
 }
 
-void DataState::ComputeDofsOffsets(std::vector<GridFunction*> &gf_array)
+std::unique_ptr<ComplexGridFunction>
+DataState::ProjectVectorFEGridFunction(std::unique_ptr<ComplexGridFunction> gf)
+{
+   if ((gf->VectorDim() == 3) && (gf->FESpace()->GetVDim() == 1))
+   {
+      int p = gf->FESpace()->GetOrder(0);
+      cout << "Switching to order " << p
+           << " discontinuous complex vector grid function..." << endl;
+      int dim = gf->FESpace()->GetMesh()->Dimension();
+      FiniteElementCollection *d_fec =
+         (p == 1 && dim == 3) ?
+         (FiniteElementCollection*)new LinearDiscont3DFECollection :
+         (FiniteElementCollection*)new L2_FECollection(p, dim, 1);
+      FiniteElementSpace *d_fespace =
+         new FiniteElementSpace(gf->FESpace()->GetMesh(), d_fec, 3);
+      ComplexGridFunction *d_gf = new ComplexGridFunction(d_fespace);
+      d_gf->MakeOwner(d_fec);
+      gf->real().ProjectVectorFieldOn(d_gf->real());
+      gf->imag().ProjectVectorFieldOn(d_gf->imag());
+      gf.reset(d_gf);
+   }
+   return gf;
+}
+
+void DataState::ProjectVectorFEGridFunction()
+{
+   if (cgrid_f)
+   {
+      internal.cgrid_f = ProjectVectorFEGridFunction(std::move(internal.cgrid_f));
+      SetComplexSolution(GetComplexSolution(), false);
+   }
+   else
+   {
+      internal.grid_f = ProjectVectorFEGridFunction(std::move(internal.grid_f));
+   }
+}
+
+void DataState::FindComplexValueRange(double &minv, double &maxv,
+                                      std::function<double(double)> scale) const
+{
+   if (!cgrid_f)
+   {
+      std::cerr << "Not a complex function" << std::endl;
+      return;
+   }
+
+   if (cmplx_sol == ComplexSolution::Magnitude)
+   {
+      // pass to the default routine
+      minv = maxv = 0.;
+      return;
+   }
+
+   if (cmplx_sol == ComplexSolution::Phase)
+   {
+      minv = -M_PI;
+      maxv = +M_PI;
+      return;
+   }
+
+   double cmplx_mag_max = 0.;
+
+   const FiniteElementSpace *fes = cgrid_f->FESpace();
+   const Mesh *msh = fes->GetMesh();
+   const int dim = msh->Dimension();
+   MFEM_VERIFY(fes->GetVectorDim() == 1, "Not a scalar function!");
+
+   if (fes->FEColl()->GetMapType(dim) == FiniteElement::VALUE)
+   {
+      for (int i = 0; i < fes->GetNDofs(); i++)
+      {
+         double rval = cgrid_f->real()(i);
+         double ival = cgrid_f->imag()(i);
+         double mag = hypot(rval, ival);
+         if (scale) { mag = scale(mag); }
+         cmplx_mag_max = max(cmplx_mag_max, mag);
+      }
+   }
+   else
+   {
+      Array<int> vdofs;
+      ElementTransformation *Tr;
+      Vector r_data, i_data;
+      Vector shape;
+      for (int z = 0; z < msh->GetNE(); z++)
+      {
+         fes->GetElementVDofs(z, vdofs);
+         cgrid_f->real().GetSubVector(vdofs, r_data);
+         cgrid_f->imag().GetSubVector(vdofs, i_data);
+
+         Tr = fes->GetElementTransformation(z);
+         const FiniteElement *fe = fes->GetFE(z);
+         const IntegrationRule &ir = fe->GetNodes();
+         const int nnp = ir.GetNPoints();
+         shape.SetSize(nnp);
+
+         for (int n = 0; n < nnp; n++)
+         {
+            const IntegrationPoint &ip = ir.IntPoint(n);
+            Tr->SetIntPoint(&ip);
+            fe->CalcPhysShape(*Tr, shape);
+            double w = shape(n);
+            double rval = r_data(n) * w;
+            double ival = i_data(n) * w;
+            double mag = hypot(rval, ival);
+            if (scale) { mag = scale(mag); }
+            cmplx_mag_max = max(cmplx_mag_max, mag);
+         }
+      }
+   }
+
+   switch (cmplx_sol)
+   {
+      case ComplexSolution::Real:
+      case ComplexSolution::Imag:
+         minv = -cmplx_mag_max;
+         maxv = +cmplx_mag_max;
+         break;
+      default:
+         std::cerr << "Unkown complex representation" << std::endl;
+         break;
+   }
+}
+
+void DataState::FindComplexValueRange(double &minv, double &maxv,
+                                      std::function<double(const Vector &)> vec2scalar,
+                                      std::function<double(double)> scale) const
+{
+   if (!cgrid_f)
+   {
+      std::cerr << "Not a complex function" << std::endl;
+      return;
+   }
+
+   if (cmplx_sol == ComplexSolution::Magnitude)
+   {
+      // pass to the default routine
+      minv = maxv = 0.;
+      return;
+   }
+
+   const int dim = mesh->Dimension();
+   const int vdim = cgrid_f->VectorDim();
+
+   // detect common symmetries
+
+   enum class Symmetry { Unknown, Norm, Phase, Linear };
+
+   Symmetry sym = Symmetry::Unknown;
+
+   {
+      double a, b;
+      Vector va(vdim), vb(vdim);
+      va = -1.;
+      vb = +1.;
+      a = vec2scalar(va);
+      b = vec2scalar(vb);
+      if (a != 0. && b != 0.)
+      {
+         if (a == b) { sym = Symmetry::Norm; }
+         else if (a == -b) { sym = Symmetry::Linear; }
+         else { sym = Symmetry::Phase; } // most likely
+      }
+   }
+
+   if (sym == Symmetry::Unknown)
+   {
+      // pass to the default routine
+      minv = maxv = 0.;
+      return;
+   }
+
+   // output phase
+   if (sym == Symmetry::Phase)
+   {
+      minv = -M_PI;
+      maxv = +M_PI;
+      return;
+   }
+
+   // complex phase
+   if (cmplx_sol == ComplexSolution::Phase)
+   {
+      Vector va(vdim), vb(vdim);
+      va = -M_PI;
+      vb = +M_PI;
+      minv = (sym == Symmetry::Norm)?(0.):(vec2scalar(va));
+      maxv = vec2scalar(vb);
+      return;
+   }
+
+   // real or imag components
+
+   function<void(const Vector &rd, const Vector &id)> kernel;
+   if (sym == Symmetry::Norm)
+   {
+      kernel = [&](const Vector &rd, const Vector &id)
+      {
+         Vector md(vdim);
+         for (int d = 0; d < vdim; d++) { md(d) = hypot(rd(d), id(d)); }
+         double scal = vec2scalar(md);
+         if (scale) { scal = scale(scal); }
+         minv = min(minv, scal);
+         maxv = max(maxv, scal);
+      };
+   }
+   else
+   {
+      kernel = [&](const Vector &rd, const Vector &id)
+      {
+         Vector md(vdim);
+         for (int d = 0; d < vdim; d++) { md(d) = hypot(rd(d), id(d)); }
+         double scal = vec2scalar(md);
+         if (scale) { scal = scale(scal); }
+         maxv = max(maxv, scal);
+      };
+   }
+
+   minv = INFINITY;
+   maxv = 0.;
+
+   const FiniteElementSpace *fes = cgrid_f->FESpace();
+
+   if (fes->FEColl()->GetMapType(dim) == FiniteElement::VALUE
+       && fes->FEColl()->GetRangeType(dim) == FiniteElement::SCALAR)
+   {
+      const int ndofs = fes->GetNDofs();
+      Array<int> vdofs;
+      Vector rd, id;
+      for (int i = 0; i < ndofs; i++)
+      {
+         vdofs.SetSize(1);
+         vdofs[0] = i;
+         fes->DofsToVDofs(vdofs);
+         cgrid_f->real().GetSubVector(vdofs, rd);
+         cgrid_f->imag().GetSubVector(vdofs, id);
+         kernel(rd, id);
+      }
+   }
+   else
+   {
+      Array<int> vdofs;
+      ElementTransformation *Tr;
+      Vector r_data, i_data;
+      Vector r_vec(vdim), i_vec(vdim);
+      DenseMatrix vshape;
+      Vector shape;
+
+      for (int z = 0; z < mesh->GetNE(); z++)
+      {
+         fes->GetElementVDofs(z, vdofs);
+         cgrid_f->real().GetSubVector(vdofs, r_data);
+         cgrid_f->imag().GetSubVector(vdofs, i_data);
+
+         Tr = fes->GetElementTransformation(z);
+         const FiniteElement *fe = fes->GetFE(z);
+         const IntegrationRule &ir = fe->GetNodes();
+         const int nnp = ir.GetNPoints();
+
+         if (fe->GetRangeType() == FiniteElement::SCALAR)
+         {
+            shape.SetSize(nnp);
+
+            for (int n = 0; n < nnp; n++)
+            {
+               const IntegrationPoint &ip = ir.IntPoint(n);
+               Tr->SetIntPoint(&ip);
+               fe->CalcPhysShape(*Tr, shape);
+               double w = shape(n);
+               for (int d = 0; d < vdim; d++)
+               {
+                  r_vec(d) = r_data(n+d*nnp) * w;
+                  i_vec(d) = i_data(n+d*nnp) * w;
+               }
+               kernel(r_vec, i_vec);
+            }
+         }
+         else
+         {
+            vshape.SetSize(nnp, vdim);
+            for (int n = 0; n < nnp; n++)
+            {
+               const IntegrationPoint &ip = ir.IntPoint(n);
+               Tr->SetIntPoint(&ip);
+               fe->CalcPhysVShape(*Tr, vshape);
+               for (int d = 0; d < vdim; d++)
+               {
+                  r_vec(d) = r_data(n) * vshape(n,d);
+                  i_vec(d) = i_data(n) * vshape(n,d);
+               }
+               kernel(r_vec, i_vec);
+            }
+         }
+      }
+   }
+
+   if (sym != Symmetry::Norm) { minv = -maxv; }
+}
+
+void DataState::ComputeDofsOffsets(std::vector<const FiniteElementSpace*>
+                                   &fespaces)
 {
-   const int nprocs = static_cast<int>(gf_array.size());
-   MFEM_VERIFY(!gf_array.empty(), "No grid functions provided for offsets");
+   const int nprocs = static_cast<int>(fespaces.size());
+   MFEM_VERIFY(!fespaces.empty(), "No grid functions provided for offsets");
 
    // only 2D meshes are supported for dofs offsets computation
-   if (gf_array[0]->FESpace()->GetMesh()->Dimension() != 2) { return; }
+   if (fespaces[0]->GetMesh()->Dimension() != 2) { return; }
 
    internal.offsets = std::make_unique<DataState::Offsets>(nprocs);
 
@@ -640,8 +1286,7 @@ void DataState::ComputeDofsOffsets(std::vector<GridFunction*> &gf_array)
    Array<int> dofs, vertices;
    for (int rank = 0, g_e = 0; rank < nprocs; rank++)
    {
-      const GridFunction *gf = gf_array[rank];
-      const FiniteElementSpace *l_fes = gf->FESpace();
+      const FiniteElementSpace *l_fes = fespaces[rank];
       Mesh *l_mesh = l_fes->GetMesh();
       // store the dofs numbers as they are fespace dependent
       auto &offset = (*offsets)[rank];
diff --git a/lib/data_state.hpp b/lib/data_state.hpp
index a5d192c8..47a0d5bc 100644
--- a/lib/data_state.hpp
+++ b/lib/data_state.hpp
@@ -17,6 +17,7 @@
 #include <memory>
 #include <vector>
 #include <utility>
+#include <functional>
 
 #include <mfem.hpp>
 
@@ -48,6 +49,19 @@ struct DataState
       MAX
    };
 
+   enum class ComplexSolution
+   {
+      NONE = -1,
+      MIN = -1,
+      //----------
+      Magnitude,
+      Phase,
+      Real,
+      Imag,
+      //----------
+      MAX
+   };
+
    // Class used for storing offsets and map of DOFs for each rank
    class Offset
    {
@@ -76,19 +90,35 @@ struct DataState
       std::unique_ptr<mfem::Mesh> mesh;
       std::unique_ptr<mfem::Mesh> mesh_quad;
       std::unique_ptr<mfem::GridFunction> grid_f;
+      std::unique_ptr<mfem::ComplexGridFunction> cgrid_f;
       std::unique_ptr<mfem::QuadratureFunction> quad_f;
       std::unique_ptr<mfem::DataCollection> data_coll;
       std::unique_ptr<Offsets> offsets;
    } internal;
 
    FieldType type {FieldType::UNKNOWN};
+   ComplexSolution cmplx_sol {ComplexSolution::NONE};
    QuadSolution quad_sol {QuadSolution::NONE};
 
    void SetGridFunctionSolution(int component = -1);
+   void SetComplexFunctionSolution(int component = -1);
    void SetQuadFunctionSolution(int component = -1);
 
+   static std::unique_ptr<mfem::GridFunction>
+   ProjectVectorFEGridFunction(std::unique_ptr<mfem::GridFunction> gf);
+
+   static std::unique_ptr<mfem::ComplexGridFunction>
+   ProjectVectorFEGridFunction(std::unique_ptr<mfem::ComplexGridFunction> gf);
+
+   void FindComplexValueRange(double &minv, double &maxv,
+                              std::function<double(double)> = {}) const;
+
+   void FindComplexValueRange(double &minv, double &maxv,
+                              std::function<double(const mfem::Vector &)>,
+                              std::function<double(double)> = {}) const;
+
    /// Compute the dofs offsets from the grid function vector
-   void ComputeDofsOffsets(std::vector<mfem::GridFunction*> &gf_array);
+   void ComputeDofsOffsets(std::vector<const mfem::FiniteElementSpace*> &fespaces);
 
 public:
    const std::unique_ptr<mfem::Vector> &sol{internal.sol};
@@ -99,6 +129,7 @@ struct DataState
    const std::unique_ptr<mfem::Mesh> &mesh{internal.mesh};
    const std::unique_ptr<mfem::Mesh> &mesh_quad{internal.mesh_quad};
    const std::unique_ptr<mfem::GridFunction> &grid_f{internal.grid_f};
+   const std::unique_ptr<mfem::ComplexGridFunction> &cgrid_f{internal.cgrid_f};
    const std::unique_ptr<mfem::QuadratureFunction> &quad_f{internal.quad_f};
    const std::unique_ptr<mfem::DataCollection> &data_coll{internal.data_coll};
    const std::unique_ptr<Offsets> &offsets{internal.offsets};
@@ -107,6 +138,7 @@ struct DataState
    bool fix_elem_orient{false};
    bool save_coloring{false};
    bool keep_attr{false};
+   double cmplx_phase{0.};
 
    DataState() = default;
    DataState(DataState &&ss) { *this = std::move(ss); }
@@ -153,6 +185,23 @@ struct DataState
    void SetGridFunction(std::vector<mfem::GridFunction*> &gf_array,
                         int num_pieces, int component = -1);
 
+   /// Set a complex grid function (plain pointer version)
+   /** Note that ownership is passed from the caller.
+       @see SetCmplxGridFunction(std::unique_ptr<mfem::ComplexGridFunction> &&, int ) */
+   void SetCmplxGridFunction(mfem::ComplexGridFunction *gf, int component = -1);
+
+   /// Set a complex grid function (unique pointer version)
+   /** Sets the complex grid function or its component (-1 means all
+       components). */
+   void SetCmplxGridFunction(std::unique_ptr<mfem::ComplexGridFunction> &&pgf,
+                             int component = -1);
+
+   /// Set a complex grid function from pieces
+   /** Serializes the pieces of a complex grid function and sets it or its
+       component (-1 means all components) */
+   void SetCmplxGridFunction(const std::vector<mfem::ComplexGridFunction*>
+                             &cgf_array, int component = -1);
+
    /// Set a quadrature function (plain pointer version)
    /** Note that ownership is passed from the caller.
        @see SetQuadFunction(std::unique_ptr<mfem::QuadFunction> &&, int ) */
@@ -192,6 +241,13 @@ struct DataState
    /// Set a (checkerboard) solution when only the mesh is given
    void SetMeshSolution();
 
+   /// Set the complex function representation producing a proxy grid function
+   void SetComplexSolution(ComplexSolution type = ComplexSolution::Magnitude,
+                           bool print = true);
+
+   /// Get the current representation of complex solution
+   inline ComplexSolution GetComplexSolution() const { return cmplx_sol; }
+
    /// Set the quadrature function representation producing a proxy grid function
    void SetQuadSolution(QuadSolution type = QuadSolution::LOR_ClosedGL);
 
@@ -204,11 +260,22 @@ struct DataState
    // Replace a given VectorFiniteElement-based grid function (e.g. from a Nedelec
    // or Raviart-Thomas space) with a discontinuous piece-wise polynomial Cartesian
    // product vector grid function of the same order.
-   static std::unique_ptr<mfem::GridFunction>
-   ProjectVectorFEGridFunction(std::unique_ptr<mfem::GridFunction> gf);
-
-   void ProjectVectorFEGridFunction()
-   { internal.grid_f = ProjectVectorFEGridFunction(std::move(internal.grid_f)); }
+   void ProjectVectorFEGridFunction();
+
+   /// Find value range of the scalar data for an intermediate representation
+   /** Returns @p minv = @p maxv = 0 if the range should be normally determined
+       from the grid function representation. */
+   void FindValueRange(double &minv, double &maxv,
+                       std::function<double(double)> scale = {}) const
+   { if (cgrid_f) { FindComplexValueRange(minv, maxv, scale); } else { minv = maxv = 0.; } }
+
+   /// Find value range of the vector data for an intermediate representation
+   /** Returns @p minv = @p maxv = 0 if the range should be normally determined
+       from the grid function representation. */
+   void FindValueRange(double &minv, double &maxv,
+                       std::function<double(const mfem::Vector &)> vec2scal,
+                       std::function<double(double)> scale = {}) const
+   { if (cgrid_f) { FindComplexValueRange(minv, maxv, vec2scal, scale); } else { minv = maxv = 0.; } }
 };
 
 #endif // GLVIS_DATA_STATE_HPP
diff --git a/lib/file_reader.cpp b/lib/file_reader.cpp
index 7e8b3a04..b4844a01 100644
--- a/lib/file_reader.cpp
+++ b/lib/file_reader.cpp
@@ -9,10 +9,11 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#include <vector>
-
 #include "file_reader.hpp"
 
+#include <vector>
+#include <general/text.hpp>
+
 using namespace std;
 using namespace mfem;
 
@@ -54,8 +55,18 @@ int FileReader::ReadSerial(FileReader::FileType ft, const char *mesh_file,
       switch (ft)
       {
          case FileType::GRID_FUNC:
-            data.SetGridFunction(new GridFunction(data.mesh.get(), *solin), component);
-            break;
+         {
+            if (CheckStreamIsComplex(*solin))
+            {
+               data.SetCmplxGridFunction(new ComplexGridFunction(data.mesh.get(), *solin),
+                                         component);
+            }
+            else
+            {
+               data.SetGridFunction(new GridFunction(data.mesh.get(), *solin), component);
+            }
+         }
+         break;
          case FileType::QUAD_FUNC:
             data.SetQuadFunction(new QuadratureFunction(data.mesh.get(), *solin),
                                  component);
@@ -125,6 +136,19 @@ int FileReader::ReadParallel(int np, FileType ft, const char *mesh_file,
    return read_err;
 }
 
+bool FileReader::CheckStreamIsComplex(std::istream &solin, bool parallel)
+{
+   string buff;
+   auto pos = solin.tellg();
+   solin >> ws;
+   getline(solin, buff);
+   solin.seekg(pos);
+   mfem::filter_dos(buff);
+   const char *header = (parallel)?("ParComplexGridFunction"):
+                        ("ComplexGridFunction");
+   return (buff == header);
+}
+
 int FileReader::ReadParMeshAndGridFunction(int np, const char *mesh_prefix,
                                            const char *sol_prefix, int component)
 {
@@ -140,6 +164,10 @@ int FileReader::ReadParMeshAndGridFunction(int np, const char *mesh_prefix,
 
    std::vector<Mesh *> mesh_array(np);
    std::vector<GridFunction *> gf_array(np);
+   std::vector<ComplexGridFunction *> cgf_array(np);
+
+   int gf_count = 0;
+   int cgf_count = 0;
 
    int read_err = 0;
    for (int p = 0; p < np; p++)
@@ -185,22 +213,59 @@ int FileReader::ReadParMeshAndGridFunction(int np, const char *mesh_prefix,
                break;
             }
 
-            gf_array[p] = new GridFunction(mesh_array[p], solfile);
+            if (CheckStreamIsComplex(solfile, true))
+            {
+               solfile >> ws;
+               solfile.ignore(3);// ignore 'Par' prefix to load as serial
+               cgf_array[p] = new ComplexGridFunction(mesh_array[p], solfile);
+               cgf_count++;
+            }
+            else
+            {
+               gf_array[p] = new GridFunction(mesh_array[p], solfile);
+               gf_count++;
+            }
          }
          else  // mesh and solution in the same file
          {
-            gf_array[p] = new GridFunction(mesh_array[p], meshfile);
+            if (CheckStreamIsComplex(meshfile))
+            {
+               cgf_array[p] = new ComplexGridFunction(mesh_array[p], meshfile);
+               cgf_count++;
+            }
+            else
+            {
+               gf_array[p] = new GridFunction(mesh_array[p], meshfile);
+               gf_count++;
+            }
          }
       }
    }
 
+   if (!read_err)
+   {
+      if ((gf_count > 0 && gf_count != np)
+          || (cgf_count > 0 && cgf_count != np))
+      {
+         cerr << "Input files contain a mixture of data types!" << endl;
+         read_err = 3;
+      }
+   }
+
    if (!read_err)
    {
       // create the combined mesh and gf
       data.SetMesh(new Mesh(mesh_array.data(), np));
       if (sol_prefix)
       {
-         data.SetGridFunction(gf_array, np, component);
+         if (cgf_array[0])
+         {
+            data.SetCmplxGridFunction(cgf_array, component);
+         }
+         else
+         {
+            data.SetGridFunction(gf_array, np, component);
+         }
       }
       else
       {
@@ -211,6 +276,7 @@ int FileReader::ReadParMeshAndGridFunction(int np, const char *mesh_prefix,
    for (int p = 0; p < np; p++)
    {
       delete gf_array[np-1-p];
+      delete cgf_array[np-1-p];
       delete mesh_array[np-1-p];
    }
 
diff --git a/lib/file_reader.hpp b/lib/file_reader.hpp
index 87eed714..ce176df0 100644
--- a/lib/file_reader.hpp
+++ b/lib/file_reader.hpp
@@ -19,6 +19,8 @@ class FileReader
    DataState &data;
    int pad_digits;
 
+   static bool CheckStreamIsComplex(std::istream &sol, bool parallel = false);
+
    int ReadParMeshAndGridFunction(int np, const char *mesh_prefix,
                                   const char *sol_prefix, int component = -1);
    int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix,
diff --git a/lib/openglvis.hpp b/lib/openglvis.hpp
index 0311d13e..1fe544e5 100644
--- a/lib/openglvis.hpp
+++ b/lib/openglvis.hpp
@@ -209,6 +209,8 @@ class VisualizationScene
    void Scale(double s);
    void Scale(double s1, double s2, double s3);
 
+   virtual void UpdateComplexPhase(double phase) { }
+
    void CenterObject();
    void CenterObject2D();
 
diff --git a/lib/script_controller.cpp b/lib/script_controller.cpp
index 6ee814bd..7a30778a 100644
--- a/lib/script_controller.cpp
+++ b/lib/script_controller.cpp
@@ -645,19 +645,19 @@ bool ScriptController::ExecuteScriptCommand()
             cout << "Script: autoscale: " << word;
             if (word == "off")
             {
-               win.vs->SetAutoscale(0);
+               win.vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::None);
             }
             else if (word == "on")
             {
-               win.vs->SetAutoscale(1);
+               win.vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::MeshAndValue);
             }
             else if (word == "value")
             {
-               win.vs->SetAutoscale(2);
+               win.vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::Value);
             }
             else if (word == "mesh")
             {
-               win.vs->SetAutoscale(3);
+               win.vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::Mesh);
             }
             else
             {
diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp
index c9207b92..ebab0451 100644
--- a/lib/stream_reader.cpp
+++ b/lib/stream_reader.cpp
@@ -13,6 +13,7 @@
 
 #include <array>
 #include <algorithm>
+#include <general/text.hpp>
 
 using namespace std;
 using namespace mfem;
@@ -102,6 +103,18 @@ bool StreamReader::SupportsDataType(const string &data_type)
    return it != commands.end();
 }
 
+bool StreamReader::CheckStreamIsComplex(std::istream &solin,
+                                        bool parallel)
+{
+   const char *header = (parallel)?("ParComplexGridFunction"):
+                        ("ComplexGridFunction");
+   solin >> ws;
+   // Returning of the characters to stream(buffer) does not work reliably,
+   // so compare only the initial character, which fortunately does not
+   // conincide with the header of FiniteElementSpace.
+   return (solin.peek() == header[0]);
+}
+
 int StreamReader::ReadStream(
    istream &is, const string &data_type)
 {
@@ -166,7 +179,14 @@ int StreamReader::ReadStream(
       case Command::VFem3D_GF_keys:
       case Command::Solution:
          data.SetMesh(new Mesh(is, 1, 0, data.fix_elem_orient));
-         data.SetGridFunction(new GridFunction(data.mesh.get(), is));
+         if (CheckStreamIsComplex(is))
+         {
+            data.SetCmplxGridFunction(new ComplexGridFunction(data.mesh.get(), is));
+         }
+         else
+         {
+            data.SetGridFunction(new GridFunction(data.mesh.get(), is));
+         }
          break;
       case Command::Quadrature:
          data.SetMesh(new Mesh(is, 1, 0, data.fix_elem_orient));
@@ -311,11 +331,13 @@ int StreamReader::ReadStreams(const StreamCollection &input_streams)
    const int nproc = input_streams.size();
    std::vector<Mesh*> mesh_array(nproc);
    std::vector<GridFunction*> gf_array(nproc);
+   std::vector<ComplexGridFunction*> cgf_array(nproc);
    std::vector<QuadratureFunction*> qf_array(nproc);
 
    std::string data_type;
 
    int gf_count = 0;
+   int cgf_count = 0;
    int qf_count = 0;
 
    for (int p = 0; p < nproc; p++)
@@ -350,8 +372,17 @@ int StreamReader::ReadStreams(const StreamCollection &input_streams)
       }
       else if (data_type != "mesh")
       {
-         gf_array[p] = new GridFunction(mesh_array[p], isock);
-         gf_count++;
+         if (CheckStreamIsComplex(isock, true))
+         {
+            isock.ignore(3);// ignore 'Par' prefix to load as serial
+            cgf_array[p] = new ComplexGridFunction(mesh_array[p], isock);
+            cgf_count++;
+         }
+         else
+         {
+            gf_array[p] = new GridFunction(mesh_array[p], isock);
+            gf_count++;
+         }
       }
 #ifdef GLVIS_DEBUG
       cout << "done." << endl;
@@ -359,6 +390,7 @@ int StreamReader::ReadStreams(const StreamCollection &input_streams)
    }
 
    if ((gf_count > 0 && gf_count != nproc)
+       || (cgf_count > 0 && cgf_count != nproc)
        || (qf_count > 0 && qf_count != nproc))
    {
       mfem_error("Input streams contain a mixture of data types!");
@@ -369,6 +401,10 @@ int StreamReader::ReadStreams(const StreamCollection &input_streams)
    {
       data.SetGridFunction(gf_array, nproc);
    }
+   else if (cgf_count > 0)
+   {
+      data.SetCmplxGridFunction(cgf_array);
+   }
    else if (qf_count > 0)
    {
       data.SetQuadFunction(qf_array);
@@ -405,6 +441,12 @@ void StreamReader::WriteStream(std::ostream &os)
       }
       data.quad_f->Save(os);
    }
+   else if (data.cgrid_f)
+   {
+      os << "solution\n";
+      data.mesh->Print(os);
+      data.cgrid_f->Save(os);
+   }
    else if (data.grid_f)
    {
       os << "solution\n";
diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp
index b62211af..b4ac669c 100644
--- a/lib/stream_reader.hpp
+++ b/lib/stream_reader.hpp
@@ -24,6 +24,8 @@ class StreamReader
 {
    DataState &data;
 
+   static bool CheckStreamIsComplex(std::istream &sol, bool parallel = false);
+
 public:
 
    StreamReader(DataState &data_): data(data_) { }
diff --git a/lib/threads.cpp b/lib/threads.cpp
index 67b14f92..1680a840 100644
--- a/lib/threads.cpp
+++ b/lib/threads.cpp
@@ -493,6 +493,19 @@ int GLVisCommand::Execute()
                   }
                   new_state.ExtrudeMeshAndSolution();
                }
+               else if (new_state.cgrid_f)
+               {
+                  auto cs = win.data_state.GetComplexSolution();
+                  if (cs != DataState::ComplexSolution::NONE)
+                  {
+                     new_state.SetComplexSolution(cs);
+                  }
+                  else
+                  {
+                     new_state.SetComplexSolution();
+                  }
+                  new_state.ExtrudeMeshAndSolution();
+               }
                break;
             default:
                cerr << "Unknown field type" << endl;
@@ -693,19 +706,19 @@ int GLVisCommand::Execute()
          cout << "Command: autoscale: " << autoscale_mode;
          if (autoscale_mode == "off")
          {
-            win.vs->SetAutoscale(0);
+            win.vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::None);
          }
          else if (autoscale_mode == "on")
          {
-            win.vs->SetAutoscale(1);
+            win.vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::MeshAndValue);
          }
          else if (autoscale_mode == "value")
          {
-            win.vs->SetAutoscale(2);
+            win.vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::Value);
          }
          else if (autoscale_mode == "mesh")
          {
-            win.vs->SetAutoscale(3);
+            win.vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::Mesh);
          }
          else
          {
diff --git a/lib/vsdata.cpp b/lib/vsdata.cpp
index c811e5ac..8f7f8e53 100644
--- a/lib/vsdata.cpp
+++ b/lib/vsdata.cpp
@@ -103,15 +103,15 @@ int VisualizationSceneScalarData::GetAutoRefineFactor()
 
 void VisualizationSceneScalarData::DoAutoscale(bool prepare)
 {
-   if (autoscale == 1)
+   if (autoscale == Autoscale::MeshAndValue)
    {
       FindNewBoxAndValueRange(prepare);
    }
-   else if (autoscale == 2)
+   else if (autoscale == Autoscale::Value)
    {
       FindNewValueRange(prepare);
    }
-   else if (autoscale == 3)
+   else if (autoscale == Autoscale::Mesh)
    {
       FindMeshBox(prepare);
    }
@@ -119,7 +119,7 @@ void VisualizationSceneScalarData::DoAutoscale(bool prepare)
 
 void VisualizationSceneScalarData::DoAutoscaleValue(bool prepare)
 {
-   if (autoscale == 1 || autoscale == 3)
+   if (autoscale == Autoscale::MeshAndValue || autoscale == Autoscale::Mesh)
    {
       FindNewBoxAndValueRange(prepare);
    }
@@ -664,11 +664,12 @@ void Key_Mod_a_Pressed(GLenum state)
 {
    if (state & KMOD_CTRL)
    {
-      static const char *autoscale_modes[] = { "off", "on", "value", "mesh" };
-      int autoscale = vsdata->GetAutoscale();
-      autoscale = (autoscale + 1)%4;
+      constexpr int nmodes = (int) VisualizationSceneScalarData::Autoscale::Max;
+      static const char *autoscale_modes[nmodes] = { "off", "on", "value", "mesh" };
+      int autoscale = (int)vsdata->GetAutoscale();
+      autoscale = (autoscale + 1) % (nmodes);
       cout << "Autoscale: " << flush;
-      vsdata->SetAutoscale(autoscale);
+      vsdata->SetAutoscale((VisualizationSceneScalarData::Autoscale) autoscale);
       cout << autoscale_modes[autoscale] << endl;
       SendExposeEvent();
    }
@@ -1344,12 +1345,16 @@ void VisualizationSceneScalarData::ToggleTexture()
    }
 }
 
-void VisualizationSceneScalarData::SetAutoscale(int _autoscale)
+void VisualizationSceneScalarData::SetAutoscale(Autoscale _autoscale,
+                                                bool update)
 {
    if (autoscale != _autoscale)
    {
       autoscale = _autoscale;
-      DoAutoscale(true);
+      if (update)
+      {
+         DoAutoscale(true);
+      }
    }
 }
 
@@ -1469,7 +1474,7 @@ void VisualizationSceneScalarData::Init()
 
    PrepareRuler();
 
-   autoscale = 1;
+   autoscale = VisualizationSceneScalarData::Autoscale::MeshAndValue;
 }
 
 VisualizationSceneScalarData::~VisualizationSceneScalarData()
diff --git a/lib/vsdata.hpp b/lib/vsdata.hpp
index f7ad8481..3c55d3e7 100644
--- a/lib/vsdata.hpp
+++ b/lib/vsdata.hpp
@@ -62,6 +62,20 @@ class VisualizationSceneScalarData : public VisualizationScene
       Max
    };
 
+   // autoscale controls the behavior when the mesh/solution are updated
+   enum class Autoscale
+   {
+      Invalid = -1,
+      Min = -1,
+      //---------
+      None, // do not change the bounding box and the value range
+      MeshAndValue, // recompute both the bounding box and the value range (default)
+      Value, // recompute only the value range
+      Mesh, // recompute only the bounding box
+      //---------
+      Max
+   };
+
 protected:
    mfem::Mesh   *mesh{}, *mesh_coarse{};
    mfem::Vector *sol{};
@@ -102,12 +116,7 @@ class VisualizationSceneScalarData : public VisualizationScene
    int ruler_on;
    double ruler_x, ruler_y, ruler_z;
 
-   // autoscale controls the behavior when the mesh/solution are updated:
-   // 0 - do not change the bounding box and the value range
-   // 1 - recompute both the bounding box and the value range (default)
-   // 2 - recompute only the value range
-   // 3 - recompute only the bounding box
-   int autoscale;
+   Autoscale autoscale;
 
    bool logscale;
 
@@ -188,10 +197,10 @@ class VisualizationSceneScalarData : public VisualizationScene
    virtual void FindMeshBox(bool prepare) { FindNewBox(prepare); }
 
    // Perform autoscaling depending on the value of 'autoscale':
-   // 0 - do nothing
-   // 1 - call FindNewBoxAndValueRange
-   // 2 - call FindNewValueRange
-   // 3 - call FindMeshBox
+   // None - do nothing
+   // MeshAndValue - call FindNewBoxAndValueRange
+   // Value - call FindNewValueRange
+   // Mesh - call FindMeshBox
    void DoAutoscale(bool prepare);
    // Similar to the above but force recomputation of the value range
    void DoAutoscaleValue(bool prepare);
@@ -224,6 +233,8 @@ class VisualizationSceneScalarData : public VisualizationScene
 
    gl3::SceneInfo GetSceneObjs() override;
 
+   void UpdateComplexPhase(double ph) override { win.UpdateComplexPhase(ph); }
+
    void ProcessUpdatedBufs(gl3::SceneInfo& scene);
 
    void glTF_ExportBox(glTF_Builder &bld,
@@ -315,8 +326,8 @@ class VisualizationSceneScalarData : public VisualizationScene
 
    void Toggle2DView();
 
-   void SetAutoscale(int _autoscale);
-   int GetAutoscale() const { return autoscale; }
+   void SetAutoscale(Autoscale _autoscale, bool update = true);
+   Autoscale GetAutoscale() const { return autoscale; }
 
    /// Shrink the set of points towards attributes centers of gravity
    void ShrinkPoints(mfem::DenseMatrix &pointmat, int i, int fn, int di);
diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp
index 1dd31287..dfbe6341 100644
--- a/lib/vssolution.cpp
+++ b/lib/vssolution.cpp
@@ -65,9 +65,16 @@ std::string VisualizationSceneSolution::GetHelpString() const
       << "| o -  (De)refine elem. (NC shading) |" << endl
       << "| O -  Switch 'o' func. (NC shading) |" << endl
       << "| p/P  Cycle through color palettes  |" << endl
-      << "| q -  Quits                         |" << endl
-      << "| Q -  Cycle quadrature data mode    |" << endl
-      << "| r -  Reset the plot to 3D view     |" << endl
+      << "| q -  Quits                         |" << endl;
+   if (win.data_state.quad_f)
+   {
+      os << "| Q -  Cycle quadrature data mode    |" << endl;
+   }
+   else if (win.data_state.cgrid_f)
+   {
+      os << "| Q -  Cycle complex data mode       |" << endl;
+   }
+   os << "| r -  Reset the plot to 3D view     |" << endl
       << "| R -  Reset the plot to 2D view     |" << endl
       << "| s -  Turn on/off unit cube scaling |" << endl
       << "| S -  Take snapshot/Record a movie  |" << endl
@@ -99,8 +106,13 @@ std::string VisualizationSceneSolution::GetHelpString() const
       << "| *,/  Scale up/down                 |" << endl
       << "| +/-  Change z-scaling              |" << endl
       << "| . -  Start/stop spinning           |" << endl
-      << "| 0/Enter - Spinning speed and dir.  |" << endl
-      << "+------------------------------------+" << endl
+      << "| 0/Enter - Spinning speed and dir.  |" << endl;
+   if (win.data_state.cgrid_f)
+   {
+      os << "| Alt+. -  Start/stop phase anim.    |" << endl
+         << "| Alt+0/Enter - Phase anim. speed    |" << endl;
+   }
+   os << "+------------------------------------+" << endl
       << "| Mouse                              |" << endl
       << "+------------------------------------+" << endl
       << "| left   btn    - Rotation           |" << endl
@@ -1150,8 +1162,14 @@ void VisualizationSceneSolution::FindNewBox(bool prepare)
 {
    FindNewBox(bb.x, bb.y, bb.z);
 
-   minv = bb.z[0];
-   maxv = bb.z[1];
+   win.data_state.FindValueRange(minv, maxv,
+   [this](double v) { return LogVal(v); });
+
+   if (minv == 0. && maxv == 0.)
+   {
+      minv = bb.z[0];
+      maxv = bb.z[1];
+   }
 
    FixValueRange();
 
@@ -1164,11 +1182,16 @@ void VisualizationSceneSolution::FindNewBox(bool prepare)
 
 void VisualizationSceneSolution::FindNewValueRange(bool prepare)
 {
-   double rx[2], ry[2], rv[2];
+   win.data_state.FindValueRange(minv, maxv,
+   [this](double v) { return LogVal(v); });
 
-   FindNewBox(rx, ry, rv);
-   minv = rv[0];
-   maxv = rv[1];
+   if (minv == 0. && maxv == 0.)
+   {
+      double rx[2], ry[2], rv[2];
+      FindNewBox(rx, ry, rv);
+      minv = rv[0];
+      maxv = rv[1];
+   }
 
    FixValueRange();
 
diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp
index 8d81247a..8e50ab73 100644
--- a/lib/vssolution3d.cpp
+++ b/lib/vssolution3d.cpp
@@ -64,9 +64,16 @@ std::string VisualizationSceneSolution3d::GetHelpString() const
       << "| M -  Toggle the mesh in the CP     |" << endl
       << "| o/O  (De)refine elem, disc shading |" << endl
       << "| p/P  Cycle through color palettes  |" << endl
-      << "| q -  Quits                         |" << endl
-      << "| Q -  Cycle quadrature data mode    |" << endl
-      << "| r -  Reset the plot to 3D view     |" << endl
+      << "| q -  Quits                         |" << endl;
+   if (win.data_state.quad_f)
+   {
+      os << "| Q -  Cycle quadrature data mode    |" << endl;
+   }
+   else if (win.data_state.cgrid_f)
+   {
+      os << "| Q -  Cycle complex data mode       |" << endl;
+   }
+   os << "| r -  Reset the plot to 3D view     |" << endl
       << "| R -  Reset the plot to 2D view     |" << endl
       << "| s -  Turn on/off unit cube scaling |" << endl
       << "| S -  Take snapshot/Record a movie  |" << endl
@@ -102,8 +109,13 @@ std::string VisualizationSceneSolution3d::GetHelpString() const
       << "| *,/  Scale up/down                 |" << endl
       << "| +/-  Change z-scaling              |" << endl
       << "| . -  Start/stop spinning           |" << endl
-      << "| 0/Enter - Spinning speed and dir.  |" << endl
-      << "+------------------------------------+" << endl
+      << "| 0/Enter - Spinning speed and dir.  |" << endl;
+   if (win.data_state.cgrid_f)
+   {
+      os << "| Alt+. -  Start/stop phase anim.    |" << endl
+         << "| Alt+0/Enter - Phase anim. speed    |" << endl;
+   }
+   os << "+------------------------------------+" << endl
       << "| Mouse                              |" << endl
       << "+------------------------------------+" << endl
       << "| left   btn    - Rotation           |" << endl
@@ -1058,19 +1070,24 @@ void VisualizationSceneSolution3d::FindNewBox(bool prepare)
 
 void VisualizationSceneSolution3d::FindNewValueRange(bool prepare)
 {
-   int map_type = (GridF) ?
-                  GridF->FESpace()->FEColl()->GetMapType(mesh->Dimension()) :
-                  FiniteElement::VALUE;
+   win.data_state.FindValueRange(minv, maxv);
 
-   if (shading < Shading::Noncomforming || map_type != (int)FiniteElement::VALUE)
-   {
-      minv = sol->Min();
-      maxv = sol->Max();
-   }
-   else
+   if (minv == 0. && maxv == 0.)
    {
-      minv = GridF->Min();
-      maxv = GridF->Max();
+      int map_type = (GridF) ?
+                     GridF->FESpace()->FEColl()->GetMapType(mesh->Dimension()) :
+                     FiniteElement::VALUE;
+
+      if (shading < Shading::Noncomforming || map_type != (int)FiniteElement::VALUE)
+      {
+         minv = sol->Min();
+         maxv = sol->Max();
+      }
+      else
+      {
+         minv = GridF->Min();
+         maxv = GridF->Max();
+      }
    }
    FixValueRange();
    UpdateValueRange(prepare);
diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp
index acc7bc45..20cc1e31 100644
--- a/lib/vsvector.cpp
+++ b/lib/vsvector.cpp
@@ -50,9 +50,16 @@ std::string VisualizationSceneVector::GetHelpString() const
       << "| o -  (De)refine elem. (NC shading) |" << endl
       << "| O -  Switch 'o' func. (NC shading) |" << endl
       << "| p/P  Cycle through color palettes  |" << endl
-      << "| q -  Quits                         |" << endl
-      << "| Q -  Cycle quadrature data mode    |" << endl
-      << "| r -  Reset the plot to 3D view     |" << endl
+      << "| q -  Quits                         |" << endl;
+   if (win.data_state.quad_f)
+   {
+      os << "| Q -  Cycle quadrature data mode    |" << endl;
+   }
+   else if (win.data_state.cgrid_f)
+   {
+      os << "| Q -  Cycle complex data mode       |" << endl;
+   }
+   os << "| r -  Reset the plot to 3D view     |" << endl
       << "| R -  Reset the plot to 2D view     |" << endl
       << "| s -  Turn on/off unit cube scaling |" << endl
       << "| S -  Take snapshot/Record a movie  |" << endl
@@ -84,8 +91,13 @@ std::string VisualizationSceneVector::GetHelpString() const
       << "| *,/  Scale up/down                 |" << endl
       << "| +/-  Change z-scaling              |" << endl
       << "| . -  Start/stop spinning           |" << endl
-      << "| 0/Enter - Spinning speed and dir.  |" << endl
-      << "+------------------------------------+" << endl
+      << "| 0/Enter - Spinning speed and dir.  |" << endl;
+   if (win.data_state.cgrid_f)
+   {
+      os << "| Alt+. -  Start/stop phase anim.    |" << endl
+         << "| Alt+0/Enter - Phase anim. speed    |" << endl;
+   }
+   os << "+------------------------------------+" << endl
       << "| Mouse                              |" << endl
       << "+------------------------------------+" << endl
       << "| left   btn    - Rotation           |" << endl
@@ -478,7 +490,7 @@ void VisualizationSceneVector::NewMeshAndSolution(
 
    VisualizationSceneSolution::NewMeshAndSolution(mesh, mesh_coarse, sol, vgf);
 
-   if (autoscale)
+   if (autoscale != VisualizationSceneScalarData::Autoscale::None)
    {
       if (Vec2Scalar == VecLength)
       {
@@ -909,6 +921,48 @@ int VisualizationSceneVector::GetFunctionAutoRefineFactor()
    return VisualizationSceneScalarData::GetFunctionAutoRefineFactor(*VecGridF);
 }
 
+void VisualizationSceneVector::FindNewBox(bool prepare)
+{
+   VisualizationSceneSolution::FindNewBox(bb.x, bb.y, bb.z);
+
+   win.data_state.FindValueRange(minv, maxv,
+   [this](const Vector &x) { return Vec2Scalar(x(0), x(1)); },
+   [this](double v) { return LogVal(v); });
+
+   if (minv == 0. && maxv == 0.)
+   {
+      minv = bb.z[0];
+      maxv = bb.z[1];
+   }
+
+   FixValueRange();
+
+   bb.z[0] = minv;
+   bb.z[1] = maxv;
+
+   SetNewScalingFromBox(); // UpdateBoundingBox minus PrepareAxes
+   UpdateValueRange(prepare);
+}
+
+void VisualizationSceneVector::FindNewValueRange(bool prepare)
+{
+   win.data_state.FindValueRange(minv, maxv,
+   [this](const Vector &x) { return Vec2Scalar(x(0), x(1)); },
+   [this](double v) { return LogVal(v); });
+
+   if (minv == 0. && maxv == 0.)
+   {
+      double rx[2], ry[2], rv[2];
+      VisualizationSceneSolution::FindNewBox(rx, ry, rv);
+      minv = rv[0];
+      maxv = rv[1];
+   }
+
+   FixValueRange();
+
+   UpdateValueRange(prepare);
+}
+
 void VisualizationSceneVector::PrepareVectorField()
 {
    int rerun;
diff --git a/lib/vsvector.hpp b/lib/vsvector.hpp
index e019fa24..6e734b54 100644
--- a/lib/vsvector.hpp
+++ b/lib/vsvector.hpp
@@ -56,6 +56,9 @@ class VisualizationSceneVector : public VisualizationSceneSolution
 
    void NewMeshAndSolution(const DataState &s) override;
 
+   void FindNewBox(bool prepare) override;
+   void FindNewValueRange(bool prepare) override;
+
    virtual ~VisualizationSceneVector();
 
    std::string GetHelpString() const override;
diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp
index 57cd4315..f3014567 100644
--- a/lib/vsvector3d.cpp
+++ b/lib/vsvector3d.cpp
@@ -59,9 +59,16 @@ std::string VisualizationSceneVector3d::GetHelpString() const
       << "| n -  Displacements step forward    |" << endl
       << "| o/O  (De)refine elem, disc shading |" << endl
       << "| p/P  Cycle through color palettes  |" << endl
-      << "| q -  Quits                         |" << endl
-      << "| Q -  Cycle quadrature data mode    |" << endl
-      << "| r -  Reset the plot to 3D view     |" << endl
+      << "| q -  Quits                         |" << endl;
+   if (win.data_state.quad_f)
+   {
+      os << "| Q -  Cycle quadrature data mode    |" << endl;
+   }
+   else if (win.data_state.cgrid_f)
+   {
+      os << "| Q -  Cycle complex data mode       |" << endl;
+   }
+   os << "| r -  Reset the plot to 3D view     |" << endl
       << "| R -  Reset the plot to 2D view     |" << endl
       << "| s -  Turn on/off unit cube scaling |" << endl
       << "| S -  Take snapshot/Record a movie  |" << endl
@@ -95,8 +102,13 @@ std::string VisualizationSceneVector3d::GetHelpString() const
       << "| *,/  Scale up/down                 |" << endl
       << "| +/-  Change z-scaling              |" << endl
       << "| . -  Start/stop spinning           |" << endl
-      << "| 0/Enter - Spinning speed and dir.  |" << endl
-      << "+------------------------------------+" << endl
+      << "| 0/Enter - Spinning speed and dir.  |" << endl;
+   if (win.data_state.cgrid_f)
+   {
+      os << "| Alt+. -  Start/stop phase anim.    |" << endl
+         << "| Alt+0/Enter - Phase anim. speed    |" << endl;
+   }
+   os << "+------------------------------------+" << endl
       << "| Mouse                              |" << endl
       << "+------------------------------------+" << endl
       << "| left   btn    - Rotation           |" << endl
@@ -274,8 +286,8 @@ void VisualizationSceneVector3d::ToggleVectorField(int i)
    PrepareVectorField();
 }
 
-static const char *scal_func_name[] =
-{"magnitude", "x-component", "y-component", "z-component"};
+const char *VisualizationSceneVector3d::scal_func_name[(int)ScalarFunction::Max]
+   = {"magnitude", "x-component", "y-component", "z-component"};
 
 void VisualizationSceneVector3d::SetScalarFunction()
 {
@@ -283,7 +295,7 @@ void VisualizationSceneVector3d::SetScalarFunction()
 
    switch (scal_func)
    {
-      case 0: // magnitude
+      case ScalarFunction::Magnitude:
          for (int i = 0; i < sol->Size(); i++)
             (*sol)(i) = sqrt((*solx)(i) * (*solx)(i) +
                              (*soly)(i) * (*soly)(i) +
@@ -304,7 +316,7 @@ void VisualizationSceneVector3d::SetScalarFunction()
             }
          }
          break;
-      case 1: // x-component
+      case ScalarFunction::Component_X:
          *sol = *solx;
          if (GridF)
             for (int i = 0; i < GridF->Size(); i++)
@@ -312,7 +324,7 @@ void VisualizationSceneVector3d::SetScalarFunction()
                (*GridF)(i) = (*VecGridF)(fes->DofToVDof(i, 0));
             }
          break;
-      case 2: // y-component
+      case ScalarFunction::Component_Y:
          *sol = *soly;
          if (GridF)
             for (int i = 0; i < GridF->Size(); i++)
@@ -320,7 +332,7 @@ void VisualizationSceneVector3d::SetScalarFunction()
                (*GridF)(i) = (*VecGridF)(fes->DofToVDof(i, 1));
             }
          break;
-      case 3: // z-component
+      case ScalarFunction::Component_Z:
          *sol = *solz;
          if (GridF)
             for (int i = 0; i < GridF->Size(); i++)
@@ -328,14 +340,17 @@ void VisualizationSceneVector3d::SetScalarFunction()
                (*GridF)(i) = (*VecGridF)(fes->DofToVDof(i, 2));
             }
          break;
+      default:
+         cerr << "Unknown scalar representation" << endl;
+         return;
    }
-   win.extra_caption = scal_func_name[scal_func];
+   win.extra_caption = scal_func_name[(int)scal_func];
 }
 
 void VisualizationSceneVector3d::ToggleScalarFunction()
 {
-   scal_func = (scal_func + 1) % 4;
-   cout << "Displaying " << scal_func_name[scal_func] << endl;
+   scal_func = (ScalarFunction)(((int)scal_func + 1) % ((int)ScalarFunction::Max));
+   cout << "Displaying " << scal_func_name[(int)scal_func] << endl;
    SetScalarFunction();
    FindNewValueRange(true);
 }
@@ -384,7 +399,7 @@ void VisualizationSceneVector3d::Init()
 
    drawdisp = 0;
    drawvector = 0;
-   scal_func = 0;
+   scal_func = ScalarFunction::Magnitude;
 
    ianim = ianimd = 0;
    ianimmax = 10;
@@ -554,6 +569,52 @@ void VisualizationSceneVector3d::NewMeshAndSolution(
    PrepareDisplacedMesh();
 }
 
+void VisualizationSceneVector3d::FindNewValueRange(bool prepare)
+{
+   function<double(const Vector &)> vec2scal;
+   switch (scal_func)
+   {
+      case ScalarFunction::Magnitude:
+         vec2scal = [](const Vector &v) { return v.Norml2(); };
+         break;
+      case ScalarFunction::Component_X:
+         vec2scal = [](const Vector &v) { return v(0); };
+         break;
+      case ScalarFunction::Component_Y:
+         vec2scal = [](const Vector &v) { return v(1); };
+         break;
+      case ScalarFunction::Component_Z:
+         vec2scal = [](const Vector &v) { return v(2); };
+         break;
+      default:
+         cerr << "Unknown scalar representation" << endl;
+         return;
+   }
+
+   win.data_state.FindValueRange(minv, maxv, vec2scal);
+
+   if (minv == 0. && maxv == 0.)
+   {
+      int map_type = (GridF) ?
+                     GridF->FESpace()->FEColl()->GetMapType(mesh->Dimension()) :
+                     FiniteElement::VALUE;
+
+      if (shading < Shading::Noncomforming || map_type != (int)FiniteElement::VALUE)
+      {
+         minv = sol->Min();
+         maxv = sol->Max();
+      }
+      else
+      {
+         minv = GridF->Min();
+         maxv = GridF->Max();
+      }
+      return;
+   }
+   FixValueRange();
+   UpdateValueRange(prepare);
+}
+
 void VisualizationSceneVector3d::PrepareFlat()
 {
    int i, j;
diff --git a/lib/vsvector3d.hpp b/lib/vsvector3d.hpp
index 9e2e83a2..346f117b 100644
--- a/lib/vsvector3d.hpp
+++ b/lib/vsvector3d.hpp
@@ -21,7 +21,22 @@ class VisualizationSceneVector3d : public VisualizationSceneSolution3d
 protected:
 
    mfem::Vector *solx, *soly, *solz;
-   int drawvector, scal_func;
+   int drawvector;
+
+   enum class ScalarFunction
+   {
+      Invalid = -1,
+      Min = -1,
+      //----------
+      Magnitude,
+      Component_X,
+      Component_Y,
+      Component_Z,
+      //----------
+      Max
+   } scal_func;
+   static const char *scal_func_name[];
+
    double mesh_volume;
    gl3::GlDrawable vector_buf;
    gl3::GlDrawable displine_buf;
@@ -47,6 +62,8 @@ class VisualizationSceneVector3d : public VisualizationSceneSolution3d
 
    void NewMeshAndSolution(const DataState &s) override;
 
+   void FindNewValueRange(bool prepare) override;
+
    virtual ~VisualizationSceneVector3d();
 
    std::string GetHelpString() const override;
diff --git a/lib/window.cpp b/lib/window.cpp
index 7ee6899e..fe5fe435 100644
--- a/lib/window.cpp
+++ b/lib/window.cpp
@@ -73,11 +73,6 @@ bool Window::GLVisInitVis(StreamCollection input_streams)
 
    locwin = this;
 
-   if (data_state.quad_f)
-   {
-      wnd->setOnKeyDown('Q', SwitchQuadSolution);
-   }
-
    double mesh_range = -1.0;
    if (field_type == DataState::FieldType::SCALAR
        || field_type == DataState::FieldType::MESH)
@@ -159,7 +154,7 @@ bool Window::GLVisInitVis(StreamCollection input_streams)
       if (mesh_range > 0.0)
       {
          vs->SetValueRange(-mesh_range, mesh_range);
-         vs->SetAutoscale(0);
+         vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::None);
       }
       if (data_state.mesh->SpaceDimension() == 2
           && field_type == DataState::FieldType::MESH)
@@ -190,12 +185,33 @@ void Window::GLVisStartVis()
    std::cout << "GLVis window closed." << std::endl;
 }
 
+void Window::SwitchComplexSolution(DataState::ComplexSolution cmplx_type)
+{
+   data_state.SetComplexSolution(cmplx_type);
+   ResetMeshAndSolution(data_state);
+}
+
 void Window::SwitchQuadSolution(DataState::QuadSolution quad_type)
 {
    data_state.SwitchQuadSolution(quad_type);
    ResetMeshAndSolution(data_state);
 }
 
+void Window::UpdateComplexPhase(double ph)
+{
+   data_state.cmplx_phase += ph;
+   data_state.cmplx_phase -= floor(data_state.cmplx_phase);
+   DataState::ComplexSolution cs = data_state.GetComplexSolution();
+   // check if magnitude is viewed, which remains the same
+   if (cs == DataState::ComplexSolution::Magnitude) { return; }
+   data_state.SetComplexSolution(cs, false);
+   // do not autoscale for animation
+   auto as = vs->GetAutoscale();
+   vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::None);
+   ResetMeshAndSolution(data_state);
+   vs->SetAutoscale(as, false);
+}
+
 bool Window::SetNewMeshAndSolution(DataState new_state)
 {
    if (new_state.mesh->SpaceDimension() == data_state.mesh->SpaceDimension() &&
@@ -227,9 +243,28 @@ void Window::ResetMeshAndSolution(DataState &ss)
    vs->NewMeshAndSolution(ss);
 }
 
-
 thread_local Window *Window::locwin = NULL;
 
+void Window::SwitchSolution()
+{
+   if (locwin->data_state.cgrid_f)
+   {
+      SwitchComplexSolution();
+   }
+   else if (locwin->data_state.quad_f)
+   {
+      SwitchQuadSolution();
+   }
+}
+
+void Window::SwitchComplexSolution()
+{
+   int ics = ((int)locwin->data_state.GetComplexSolution()+1)
+             % ((int)DataState::ComplexSolution::MAX);
+   locwin->SwitchComplexSolution((DataState::ComplexSolution)ics);
+   SendExposeEvent();
+}
+
 void Window::SwitchQuadSolution()
 {
    int iqs = ((int)locwin->data_state.GetQuadSolution()+1)
diff --git a/lib/window.hpp b/lib/window.hpp
index eeb90ce1..3391a2aa 100644
--- a/lib/window.hpp
+++ b/lib/window.hpp
@@ -61,9 +61,15 @@ struct Window
    bool GLVisInitVis(StreamCollection input_streams);
    void GLVisStartVis();
 
+   /// Switch the complex function representation and update the visualization
+   void SwitchComplexSolution(DataState::ComplexSolution type);
+
    /// Switch the quadrature function representation and update the visualization
    void SwitchQuadSolution(DataState::QuadSolution type);
 
+   /// Update complex phase of the solution and update the visualization
+   void UpdateComplexPhase(double ph);
+
    /// Sets a new mesh and solution from another DataState object, and
    /// updates the VisualizationScene with the new data.
    ///
@@ -73,10 +79,16 @@ struct Window
    /// updated.
    bool SetNewMeshAndSolution(DataState new_state);
 
+   /// Switch representation of the solution
+   static void SwitchSolution();
+
 private:
    /// Thread-local singleton for key handlers
    static thread_local Window *locwin;
 
+   /// Switch representation of the complex function
+   static void SwitchComplexSolution();
+
    /// Switch representation of the quadrature function
    static void SwitchQuadSolution();
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 1e9cf3fb..339f6f7e 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -33,6 +33,8 @@ set(stream_tests
     "ex5"
     "ex6-dofs-numbering"
     "ex9"
+    "ex22-scalar"
+    "ex22-vector"
     "klein-bottle"
     "laghos"
     "mesh-explorer"
diff --git a/tests/data b/tests/data
index bd2a11d9..f1d01ff0 160000
--- a/tests/data
+++ b/tests/data
@@ -1 +1 @@
-Subproject commit bd2a11d9eabf7dcf2d83cffd906b313e56ab7ef0
+Subproject commit f1d01ff03a5e1d2149e15adef506a803100b45ec