glvis.cpp

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443271. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the GLVis visualization tool and library. For more
// information and source code availability see http://glvis.org.
//
// GLVis is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.


// GLVis - an OpenGL visualization server based on the MFEM library


#include <limits>
#include <iostream>
#include <fstream>
#include <string>
#include <cstdio>
#include <cstring>
#include <ctime>
#include <csignal>

#include <X11/keysym.h>
#include <unistd.h>

#include "mfem.hpp"
#include "lib/visual.hpp"

using namespace std;
using namespace mfem;

const char *string_none    = "(none)";
const char *string_default = "(default)";

// Global variables for command line arguments
const char *mesh_file       = string_none;
const char *sol_file        = string_none;
const char *vec_sol_file    = string_none;
const char *gfunc_file      = string_none;
const char *arg_keys        = string_none;
int         np              = 0;
int         pad_digits      = 6;
int         gf_component    = -1;
bool        fix_elem_orient = false;
bool        save_coloring   = false;
bool        keep_attr       = false;
int         window_x        = 0; // not a command line option
int         window_y        = 0; // not a command line option
int         window_w        = 400;
int         window_h        = 350;
const char *window_title    = string_default;
const char *c_plot_caption  = string_none;
string      plot_caption;
string      extra_caption;

// Global variables
int input = 1;
Mesh *mesh = NULL;
Vector sol, solu, solv, solw, normals;
GridFunction *grid_f = NULL;
int is_gf = 0;
string keys;
VisualizationSceneScalarData *vs = NULL;

GeometryRefiner GLVisGeometryRefiner;

const char *window_titles[] = { "GLVis [scalar data]",
                                "GLVis [vector data]", "GLVis [mesh]"
                              };
istream *script = NULL;
int scr_running = 0;
int scr_level = 0;
Vector *init_nodes = NULL;
double scr_min_val, scr_max_val;

Array<istream *> input_streams;

extern char **environ;

void Set_Palette(int); // defined in palettes.hpp

void PrintSampleUsage(ostream &out);

// read the mesh and the solution from a file
void ReadSerial();

// choose grid function component and set the input flag
void SetGridFunction();

// set a (checkerboard) solution when only the mesh is given
void SetMeshSolution(Mesh *mesh, GridFunction *&grid_f, bool save_coloring);

// read the mesh and the solution from multiple files
void ReadParallel();

int ReadParMeshAndGridFunction(int np, const char *mesh_prefix,
                               const char *sol_prefix, Mesh **mesh_p,
                               GridFunction **sol_p, int keep_attr);

int ReadInputStreams();

void CloseInputStreams(bool);

GridFunction *ProjectVectorFEGridFunction(GridFunction*);

void Extrude1DMeshAndSolution(Mesh **mesh_p, GridFunction **grid_f_p,
                              Vector *sol);

// Read the content of an input stream (e.g. from socket/file)
int ReadStream(istream &is, const string &data_type)
{
   // 0 - scalar data, 1 - vector data, 2 - mesh only, (-1) - unknown
   int field_type = 0;

   delete mesh; mesh = NULL;
   delete grid_f; grid_f = NULL;
   keys.clear();
   if (data_type == "fem2d_data")
   {
      mesh = new Mesh(is, 0, 0, fix_elem_orient);
      sol.Load(is, mesh->GetNV());
   }
   else if (data_type == "vfem2d_data" || data_type == "vfem2d_data_keys")
   {
      field_type = 1;
      mesh = new Mesh(is, 0, 0, fix_elem_orient);
      solu.Load(is, mesh->GetNV());
      solv.Load(is, mesh->GetNV());
      if (data_type == "vfem2d_data_keys")
      {
         is >> keys;
      }
   }
   else if (data_type == "fem3d_data")
   {
      mesh = new Mesh(is, 0, 0, fix_elem_orient);
      sol.Load(is, mesh->GetNV());
   }
   else if (data_type == "vfem3d_data" || data_type == "vfem3d_data_keys")
   {
      field_type = 1;
      mesh = new Mesh(is, 0, 0, fix_elem_orient);
      solu.Load(is, mesh->GetNV());
      solv.Load(is, mesh->GetNV());
      solw.Load(is, mesh->GetNV());
      if (data_type == "vfem3d_data_keys")
      {
         is >> keys;
      }
   }
   else if (data_type == "fem2d_gf_data" || data_type == "fem2d_gf_data_keys")
   {
      mesh = new Mesh(is, 1, 0, fix_elem_orient);
      grid_f = new GridFunction(mesh, is);
      if (data_type == "fem2d_gf_data_keys")
      {
         is >> keys;
      }
   }
   else if (data_type == "vfem2d_gf_data" || data_type == "vfem2d_gf_data_keys")
   {
      field_type = 1;
      mesh = new Mesh(is, 1, 0, fix_elem_orient);
      grid_f = new GridFunction(mesh, is);
      if (data_type == "vfem2d_gf_data_keys")
      {
         is >> keys;
      }
   }
   else if (data_type == "fem3d_gf_data" || data_type == "fem3d_gf_data_keys")
   {
      mesh = new Mesh(is, 1, 0, fix_elem_orient);
      grid_f = new GridFunction(mesh, is);
      if (data_type == "fem3d_gf_data_keys")
      {
         is >> keys;
      }
   }
   else if (data_type == "vfem3d_gf_data" || data_type == "vfem3d_gf_data_keys")
   {
      field_type = 1;
      mesh = new Mesh(is, 1, 0, fix_elem_orient);
      grid_f = new GridFunction(mesh, is);
      if (data_type == "vfem3d_gf_data_keys")
      {
         is >> keys;
      }
   }
   else if (data_type == "solution")
   {
      mesh = new Mesh(is, 1, 0, fix_elem_orient);
      grid_f = new GridFunction(mesh, is);
      field_type = (grid_f->VectorDim() == 1) ? 0 : 1;
   }
   else if (data_type == "mesh")
   {
      mesh = new Mesh(is, 1, 0, fix_elem_orient);
      SetMeshSolution(mesh, grid_f, save_coloring);
      field_type = 2;
   }
   else if (data_type == "raw_scalar_2d")
   {
      Array<Array<double> *> vertices;
      Array<Array<int> *> elements;
      Array<int> elem_types;
      string ident;
      int num_patches, num_vert, num_elem, n;
      is >> ws >> ident; // 'patches'
      is >> num_patches;
      // cout << ident << ' ' << num_patches << endl;
      vertices.SetSize(num_patches);
      vertices = NULL;
      elements.SetSize(num_patches);
      elements = NULL;
      elem_types.SetSize(num_patches);
      elem_types = 0;
      int tot_num_vert = 0;
      int tot_num_elem = 0;
      int mesh_type = 0;
      for (int i = 0; i < num_patches; i++)
      {
         is >> ws >> ident; // 'vertices'
         is >> num_vert;
         // cout << '\n' << ident << ' ' << num_vert << endl;
         // read vertices in the format: x y z nx ny nz
         vertices[i] = new Array<double>(6*num_vert);
         Array<double> &verts = *vertices[i];
         for (int j = 0; j < verts.Size(); j++)
         {
            is >> verts[j];
         }

         is >> ws >> ident; // 'triangles' or 'quads'
         if (ident == "triangles")
         {
            n = 3, mesh_type |= 1;
         }
         else
         {
            n = 4, mesh_type |= 2;
         }
         elem_types[i] = n;
         is >> num_elem;
         // cout << ident << ' ' << num_elem << endl;
         elements[i] = new Array<int>(n*num_elem);
         Array<int> &elems = *elements[i];
         for (int j = 0; j < elems.Size(); j++)
         {
            is >> elems[j];
            elems[j] += tot_num_vert;
         }
         tot_num_vert += num_vert;
         tot_num_elem += num_elem;
      }

      mesh = new Mesh(2, tot_num_vert, tot_num_elem, 0);
      sol.SetSize(tot_num_vert);
      normals.SetSize(3*tot_num_vert);

      int v_off = 0;
      for (int i = 0; i < num_patches; i++)
      {
         Array<double> &verts = *vertices[i];
         num_vert = verts.Size()/6;
         for (int j = 0; j < num_vert; j++)
         {
            mesh->AddVertex(&verts[6*j]);
            sol(v_off) = verts[6*j+2];
            normals(3*v_off+0) = verts[6*j+3];
            normals(3*v_off+1) = verts[6*j+4];
            normals(3*v_off+2) = verts[6*j+5];
            v_off++;
         }

         n = elem_types[i];
         Array<int> &elems = *elements[i];
         num_elem = elems.Size()/n;
         // int attr = 1;
         int attr = i + 1;
         if (n == 3)
            for (int j = 0; j < num_elem; j++)
            {
               mesh->AddTriangle(&elems[3*j], attr);
            }
         else
            for (int j = 0; j < num_elem; j++)
            {
               mesh->AddQuad(&elems[4*j], attr);
            }
      }

      if (mesh_type == 1)
      {
         mesh->FinalizeTriMesh(1, 0, fix_elem_orient);
      }
      else if (mesh_type == 2)
      {
         mesh->FinalizeQuadMesh(1, 0, fix_elem_orient);
      }
      else
      {
         mfem_error("Input data contains mixture of triangles and quads!");
      }

      mesh->GenerateBoundaryElements();

      for (int i = num_patches; i > 0; )
      {
         i--;
         delete elements[i];
         delete vertices[i];
      }

      field_type = 0;
   }
   else
   {
      field_type = -1;
      cerr << "Unknown data format" << endl;
      cerr << data_type << endl;
   }

   if (field_type >= 0 && field_type <= 2)
   {
      if (grid_f)
      {
         Extrude1DMeshAndSolution(&mesh, &grid_f, NULL);
      }
      else
      {
         Extrude1DMeshAndSolution(&mesh, NULL, &sol);
      }
   }

   return field_type;
}

int InitVis(int t)
{
   const char *win_title =
      (window_title == string_default) ? window_titles[t] : window_title;

   return InitVisualization(win_title, window_x, window_y,
                            window_w, window_h);
}

// Visualize the data in the global variables mesh, sol/grid_f, etc
void StartVisualization(int field_type)
{
   if (field_type < 0 || field_type > 2)
   {
      return;
   }

   if (InitVis(field_type))
   {
      cerr << "Initializing the visualization failed." << endl;
      return;
   }

   communication_thread *comm_thread = NULL;

   if (input_streams.Size() > 0)
   {
      auxModKeyFunc(XK_space, ThreadsPauseFunc);
      glvis_command = new GLVisCommand(&vs, &mesh, &grid_f, &sol, &keep_attr,
                                       &fix_elem_orient);
      comm_thread = new communication_thread(input_streams);
   }

   double mesh_range = -1.0;
   if (field_type == 0 || field_type == 2)
   {
      if (grid_f)
      {
         grid_f->GetNodalValues(sol);
      }
      if (mesh->SpaceDimension() == 2)
      {
         VisualizationSceneSolution *vss;
         if (field_type == 2)
         {
            Set_Palette(4);
         }
         if (normals.Size() > 0)
         {
            vs = vss = new VisualizationSceneSolution(*mesh, sol, &normals);
         }
         else
         {
            vs = vss = new VisualizationSceneSolution(*mesh, sol);
         }
         if (grid_f)
         {
            vss->SetGridFunction(*grid_f);
         }
         if (field_type == 2)
         {
            vs->OrthogonalProjection = 1;
            vs->light = 0;
            vs->Zoom(1.8);
         }
      }
      else if (mesh->SpaceDimension() == 3)
      {
         VisualizationSceneSolution3d *vss;
         vs = vss = new VisualizationSceneSolution3d(*mesh, sol);
         if (grid_f)
         {
            vss->SetGridFunction(grid_f);
         }
         if (field_type == 2)
         {
            if (mesh->Dimension() == 3)
            {
               Set_Palette(4);
               // Set_Palette(11);
               // Set_Material_And_Light(4,3);
            }
            else
            {
               Set_Palette(4);
            }
            vss->ToggleDrawAxes();
            vss->ToggleDrawMesh();
         }
      }
      if (field_type == 2)
      {
         if (grid_f)
         {
            mesh_range = grid_f->Max() + 1.0;
         }
         else
         {
            mesh_range = sol.Max() + 1.0;
         }
      }
   }
   else if (field_type == 1)
   {
      if (mesh->SpaceDimension() == 2)
      {
         if (grid_f)
         {
            vs = new VisualizationSceneVector(*grid_f);
         }
         else
         {
            vs = new VisualizationSceneVector(*mesh, solu, solv);
         }
      }
      else if (mesh->SpaceDimension() == 3)
      {
         if (grid_f)
         {
            grid_f = ProjectVectorFEGridFunction(grid_f);
            vs = new VisualizationSceneVector3d(*grid_f);
         }
         else
         {
            vs = new VisualizationSceneVector3d(*mesh, solu, solv, solw);
         }
      }
   }

   if (vs)
   {
      // increase the refinement factors if visualizing a GridFunction
      if (grid_f)
      {
         vs->AutoRefine();
         vs->SetShading(2, true);
      }
      if (mesh_range > 0.0)
      {
         vs->SetValueRange(-mesh_range, mesh_range);
         vs->SetAutoscale(0);
      }
      if (mesh->SpaceDimension() == 2 && field_type == 2)
      {
         SetVisualizationScene(vs, 2, keys.c_str());
      }
      else
      {
         SetVisualizationScene(vs, 3, keys.c_str());
      }
   }

   KillVisualization(); // deletes vs
   vs = NULL;
   if (input_streams.Size() > 0)
   {
      glvis_command->Terminate();
      delete comm_thread;
      delete glvis_command;
      glvis_command = NULL;
   }
   delete grid_f; grid_f = NULL;
   delete mesh; mesh = NULL;
   cout << "GLVis window closed." << endl;
}

int ScriptReadSolution(istream &scr, Mesh **mp, GridFunction **sp)
{
   string mword,sword;

   cout << "Script: solution: " << flush;
   // read the mesh
   scr >> ws >> mword; // mesh filename (can't contain spaces)
   cout << "mesh: " << mword << "; " << flush;
   named_ifgzstream imesh(mword.c_str());
   if (!imesh)
   {
      cout << "Can not open mesh file: " << mword << endl;
      return 1;
   }
   *mp = new Mesh(imesh, 1, 0, fix_elem_orient);

   // read the solution (GridFunction)
   scr >> ws >> sword;
   if (sword == mword) // mesh and solution in the same file
   {
      cout << "solution: " << mword << endl;
      *sp = new GridFunction(*mp, imesh);
   }
   else
   {
      cout << "solution: " << sword << endl;
      ifgzstream isol(sword.c_str());
      if (!isol)
      {
         cout << "Can not open solution file: " << sword << endl;
         delete *mp; *mp = NULL;
         return 2;
      }
      *sp = new GridFunction(*mp, isol);
   }

   Extrude1DMeshAndSolution(mp, sp, NULL);

   return 0;
}

int ScriptReadParSolution(istream &scr, Mesh **mp, GridFunction **sp)
{
   int np, keep_attr, err;
   string mesh_prefix, sol_prefix;

   cout << "Script: psolution: " << flush;
   // read number of processors
   scr >> np;
   cout << "# processors: " << np << "; " << flush;
   // read the mesh prefix
   scr >> ws >> mesh_prefix; // mesh prefix (can't contain spaces)
   cout << "mesh prefix: " << mesh_prefix << "; " << flush;
   scr >> ws >> keep_attr;
   if (keep_attr)
   {
      cout << "(real attributes); " << flush;
   }
   else
   {
      cout << "(processor attributes); " << flush;
   }
   // read the solution prefix
   scr >> ws >> sol_prefix;
   cout << "solution prefix: " << sol_prefix << endl;

   err = ReadParMeshAndGridFunction(np, mesh_prefix.c_str(),
                                    sol_prefix.c_str(), mp, sp, keep_attr);
   if (!err)
   {
      Extrude1DMeshAndSolution(mp, sp, NULL);
   }
   return err;
}

int ScriptReadDisplMesh(istream &scr, Mesh **mp, GridFunction **sp)
{
   Mesh *m;
   string word;

   cout << "Script: mesh: " << flush;
   scr >> ws >> word;
   {
      named_ifgzstream imesh(word.c_str());
      if (!imesh)
      {
         cout << "Can not open mesh file: " << word << endl;
         return 1;
      }
      cout << word << endl;
      m = new Mesh(imesh, 1, 0, fix_elem_orient);
   }
   Extrude1DMeshAndSolution(&m, NULL, NULL);
   if (init_nodes == NULL)
   {
      init_nodes = new Vector;
      m->GetNodes(*init_nodes);
      delete m;
      *mp = NULL;
      *sp = NULL;
   }
   else
   {
      FiniteElementCollection  *vfec = NULL;
      FiniteElementSpace *vfes;
      GridFunction *g;
      vfes = (FiniteElementSpace *)m->GetNodalFESpace();
      if (vfes == NULL)
      {
         vfec = new LinearFECollection;
         vfes = new FiniteElementSpace(m, vfec, m->SpaceDimension());
      }

      g = new GridFunction(vfes);
      if (vfec)
      {
         g->MakeOwner(vfec);
      }
      m->GetNodes(*g);
      if (g->Size() == init_nodes->Size())
      {
         subtract(*init_nodes, *g, *g);
      }
      else
      {
         cout << "Script: incompatible meshes!" << endl;
         *g = 0.0;
      }
      *mp = m;
      *sp = g;
   }

   return 0;
}

void ExecuteScriptCommand()
{
   if (!script)
   {
      cout << "No script stream defined! (Bug?)" << endl;
      return;
   }

   istream &scr = *script;
   string word;
   int done_one_command = 0;
   while (!done_one_command)
   {
      scr >> ws;
      if (!scr.good())
      {
         cout << "End of script." << endl;
         scr_level = 0;
         return;
      }
      if (scr.peek() == '#')
      {
         getline(scr, word);
         continue;
      }
      scr >> word;
      if (word == "{")
      {
         scr_level++;
      }
      else if (word == "}")
      {
         scr_level--;
         if (scr_level < 0)
         {
            scr_level = 0;
         }
      }
      else if (word == "solution" || word == "mesh" || word == "psolution")
      {
         Mesh *new_m = NULL;
         GridFunction *new_g;

         if (word == "solution")
         {
            if (ScriptReadSolution(scr, &new_m, &new_g))
            {
               done_one_command = 1;
               continue;
            }
         }
         else if (word == "mesh")
         {
            if (ScriptReadDisplMesh(scr, &new_m, &new_g))
            {
               done_one_command = 1;
               continue;
            }
            if (new_m == NULL)
            {
               cout << "Script: unexpected 'mesh' command!" << endl;
               done_one_command = 1;
               continue;
            }
         }
         else if (word == "psolution")
         {
            if (ScriptReadParSolution(scr, &new_m, &new_g))
            {
               done_one_command = 1;
               continue;
            }
         }

         if (new_m->SpaceDimension() == mesh->SpaceDimension() &&
             new_g->VectorDim() == grid_f->VectorDim())
         {
            if (new_m->SpaceDimension() == 2)
            {
               if (new_g->VectorDim() == 1)
               {
                  VisualizationSceneSolution *vss =
                     dynamic_cast<VisualizationSceneSolution *>(vs);
                  new_g->GetNodalValues(sol);
                  vss->NewMeshAndSolution(new_m, &sol, new_g);
               }
               else
               {
                  VisualizationSceneVector *vsv =
                     dynamic_cast<VisualizationSceneVector *>(vs);
                  vsv->NewMeshAndSolution(*new_g);
               }
            }
            else
            {
               if (new_g->VectorDim() == 1)
               {
                  VisualizationSceneSolution3d *vss =
                     dynamic_cast<VisualizationSceneSolution3d *>(vs);
                  new_g->GetNodalValues(sol);
                  vss->NewMeshAndSolution(new_m, &sol, new_g);
               }
               else
               {
                  new_g = ProjectVectorFEGridFunction(new_g);
                  VisualizationSceneVector3d *vss =
                     dynamic_cast<VisualizationSceneVector3d *>(vs);
                  vss->NewMeshAndSolution(new_m, new_g);
               }
            }
            delete grid_f; grid_f = new_g;
            delete mesh; mesh = new_m;

            vs->Draw();
         }
         else
         {
            cout << "Different type of mesh / solution." << endl;
            delete new_g;
            delete new_m;
         }
      }
      else if (word == "screenshot")
      {
         scr >> ws >> word;

         cout << "Script: screenshot: " << flush;

         if (Screenshot(word.c_str(), true))
         {
            cout << "Screenshot(" << word << ") failed." << endl;
            done_one_command = 1;
            continue;
         }
         cout << "-> " << word << endl;

         if (scr_min_val > vs->GetMinV())
         {
            scr_min_val = vs->GetMinV();
         }
         if (scr_max_val < vs->GetMaxV())
         {
            scr_max_val = vs->GetMaxV();
         }
      }
      else if (word == "viewcenter")
      {
         scr >> vs->ViewCenterX >> vs->ViewCenterY;
         cout << "Script: viewcenter: "
              << vs->ViewCenterX << ' ' << vs->ViewCenterY << endl;
         MyExpose();
      }
      else if (word ==  "perspective")
      {
         scr >> ws >> word;
         cout << "Script: perspective: " << word;
         if (word == "off")
         {
            vs->OrthogonalProjection = 1;
         }
         else if (word == "on")
         {
            vs->OrthogonalProjection = 0;
         }
         else
         {
            cout << '?';
         }
         cout << endl;
         MyExpose();
      }
      else if (word ==  "light")
      {
         scr >> ws >> word;
         cout << "Script: light: " << word;
         if (word == "off")
         {
            vs->light = 0;
         }
         else if (word == "on")
         {
            vs->light = 1;
         }
         else
         {
            cout << '?';
         }
         cout << endl;
         MyExpose();
      }
      else if (word == "view")
      {
         double theta, phi;
         scr >> theta >> phi;
         cout << "Script: view: " << theta << ' ' << phi << endl;
         vs->SetView(theta, phi);
         MyExpose();
      }
      else if (word == "zoom")
      {
         double factor;
         scr >> factor;
         cout << "Script: zoom: " << factor << endl;
         vs->Zoom(factor);
         MyExpose();
      }
      else if (word == "shading")
      {
         scr >> ws >> word;
         cout << "Script: shading: " << flush;
         int s = -1;
         if (word == "flat")
         {
            s = 0;
         }
         else if (word == "smooth")
         {
            s = 1;
         }
         else if (word == "cool")
         {
            s = 2;
         }
         if (s != -1)
         {
            vs->SetShading(s, false);
            cout << word << endl;
            MyExpose();
         }
         else
         {
            cout << word << " ?" << endl;
         }
      }
      else if (word == "subdivisions")
      {
         int t, b;
         scr >> t >> b;
         cout << "Script: subdivisions: " << flush;
         vs->SetRefineFactors(t, b);
         cout << t << ' ' << b << endl;
         MyExpose();
      }
      else if (word == "valuerange")
      {
         double min, max;
         scr >> min >> max;
         cout << "Script: valuerange: " << flush;
         vs->SetValueRange(min, max);
         cout << min << ' ' << max << endl;
         MyExpose();
      }
      else if (word == "autoscale")
      {
         scr >> ws >> word;
         cout << "Script: autoscale: " << word;
         if (word == "off")
         {
            vs->SetAutoscale(0);
         }
         else if (word == "on")
         {
            vs->SetAutoscale(1);
         }
         else if (word == "value")
         {
            vs->SetAutoscale(2);
         }
         else if (word == "mesh")
         {
            vs->SetAutoscale(3);
         }
         else
         {
            cout << '?';
         }
         cout << endl;
      }
      else if (word == "window")
      {
         scr >> window_x >> window_y >> window_w >> window_h;
         cout << "Script: window: " << window_x << ' ' << window_y
              << ' ' << window_w << ' ' << window_h << endl;
         MoveResizeWindow(window_x, window_y, window_w, window_h);
         MyExpose();
      }
      else if (word == "keys")
      {
         scr >> keys;
         cout << "Script: keys: '" << keys << "'" << endl;
         // SendKeySequence(keys.c_str());
         CallKeySequence(keys.c_str());
         MyExpose();
      }
      else if (word == "palette")
      {
         int pal;
         scr >> pal;
         cout << "Script: palette: " << pal << endl;
         Set_Palette(pal-1);
         if (!GetUseTexture())
         {
            vs->EventUpdateColors();
         }
         MyExpose();
      }
      else if (word == "toggle_attributes")
      {
         Array<int> attr_list;
         cout << "Script: toggle_attributes:";
         for (scr >> ws; scr.peek() != ';'; scr >> ws)
         {
            attr_list.Append(0);
            scr >> attr_list.Last();
            if (attr_list.Size() <= 256)
            {
               cout << ' ' << attr_list.Last();
            }
            else if (attr_list.Size() == 257)
            {
               cout << " ... " << flush;
            }
         }
         scr.get(); // read the end symbol: ';'
         cout << endl;
         vs->ToggleAttributes(attr_list);
         MyExpose();
      }
      else if (word == "rotmat")
      {
         cout << "Script: rotmat:";
         for (int i = 0; i < 16; i++)
         {
            scr >> vs->rotmat[i];
            cout << ' ' << vs->rotmat[i];
         }
         cout << endl;
         MyExpose();
      }
      else if (word == "camera")
      {
         double cam[9];
         cout << "Script: camera:";
         for (int i = 0; i < 9; i++)
         {
            scr >> cam[i];
            cout << ' ' << cam[i];
         }
         cout << endl;
         vs->cam.Set(cam);
         MyExpose();
      }
      else if (word == "scale")
      {
         double scale;
         cout << "Script: scale:";
         scr >> scale;
         cout << ' ' << scale;
         cout << endl;
         vs->Scale(scale);
         MyExpose();
      }
      else if (word == "translate")
      {
         double x, y, z;
         cout << "Script: translate:";
         scr >> x >> y >> z;
         cout << ' ' << x << ' ' << y << ' ' << z;
         cout << endl;
         vs->Translate(x, y, z);
         MyExpose();
      }
      else if (word == "plot_caption")
      {
         char delim;
         scr >> ws >> delim;
         getline(scr, plot_caption, delim);
         vs->UpdateCaption(); // turn on or off the caption
         MyExpose();
      }
      else
      {
         cout << "Unknown command in script: " << word << endl;
      }

      done_one_command = 1;
   }
}

void ScriptControl();

void ScriptIdleFunc()
{
   ExecuteScriptCommand();
   if (scr_level == 0)
   {
      ScriptControl();
   }
}

void ScriptControl()
{
   if (scr_running)
   {
      scr_running = 0;
      RemoveIdleFunc(ScriptIdleFunc);
   }
   else
   {
      scr_running = 1;
      AddIdleFunc(ScriptIdleFunc);
   }
}

void PlayScript(istream &scr)
{
   string word;

   scr_min_val = numeric_limits<double>::infinity();
   scr_max_val = -scr_min_val;

   // read initializing commands
   while (1)
   {
      scr >> ws;
      if (!scr.good())
      {
         cout << "Error in script" << endl;
         return;
      }
      if (scr.peek() == '#')
      {
         getline(scr, word);
         continue;
      }
      scr >> word;
      if (word == "window")
      {
         scr >> window_x >> window_y >> window_w >> window_h;
      }
      else if (word == "solution")
      {
         if (ScriptReadSolution(scr, &mesh, &grid_f))
         {
            return;
         }

         // start the visualization
         break;
      }
      else if (word == "psolution")
      {
         if (ScriptReadParSolution(scr, &mesh, &grid_f))
         {
            return;
         }

         // start the visualization
         break;
      }
      else if (word == "mesh")
      {
         if (ScriptReadDisplMesh(scr, &mesh, &grid_f))
         {
            return;
         }
         if (mesh)
         {
            break;
         }
      }
      else
      {
         cout << "Unknown command in script: " << word << endl;
      }
   }

   scr_level = scr_running = 0;
   auxKeyFunc(XK_space, ScriptControl);
   script = &scr;
   keys.clear();

   StartVisualization((grid_f->VectorDim() == 1) ? 0 : 1);

   delete init_nodes; init_nodes = NULL;

   cout << "Script: min_val = " << scr_min_val
        << ", max_val = " << scr_max_val << endl;

   script = NULL;
}


int main (int argc, char *argv[])
{
   // variables for command line arguments
   bool        multi_session = true;   // not added as option
   bool        mac           = false;
   const char *stream_file   = string_none;
   const char *script_file   = string_none;
   const char *font_name     = string_default;
   int         portnum       = 19916;
   bool        secure        = socketstream::secure_default;
   int         multisample   = GetMultisample();
   double      line_width    = Get_LineWidth();
   double      ms_line_width = Get_MS_LineWidth();
   int         geom_ref_type = Quadrature1D::ClosedUniform;

   OptionsParser args(argc, argv);

   args.AddOption(&mesh_file, "-m", "--mesh",
                  "Mesh file to visualize.");
   args.AddOption(&gfunc_file, "-g", "--grid-function",
                  "Solution (GridFunction) file to visualize.");
   args.AddOption(&gf_component, "-gc", "--grid-function-component",
                  "Select a grid function component, [0-<num-comp>) or"
                  " -1 for all.");
   args.AddOption(&sol_file, "-s", "--scalar-solution",
                  "Scalar solution (vertex values) file to visualize.");
   args.AddOption(&vec_sol_file, "-v", "--vector-solution",
                  "Vector solution (vertex values) file to visualize.");
   args.AddOption(&np, "-np", "--num-proc",
                  "Load mesh/solution from multiple processors.");
   args.AddOption(&pad_digits, "-d", "--pad-digits",
                  "Number of digits used for processor ranks in file names.");
   args.AddOption(&script_file, "-run", "--run-script",
                  "Run a GLVis script file.");
   args.AddOption(&arg_keys, "-k", "--keys",
                  "Execute key shortcut commands in the GLVis window.");
   args.AddOption(&fix_elem_orient, "-fo", "--fix-orientations",
                  "-no-fo", "--dont-fix-orientations",
                  "Attempt to fix the orientations of inverted elements.");
   args.AddOption(&keep_attr, "-a", "--real-attributes",
                  "-ap", "--processor-attributes",
                  "When opening a parallel mesh, use the real mesh attributes"
                  " or replace them with the processor rank.");
   args.AddOption(&geom_ref_type, "-grt", "--geometry-refiner-type",
                  "Set of points to use when refining geometry:"
                  " 3 = uniform, 1 = Gauss-Lobatto, (see mfem::Quadrature1D).");
   args.AddOption(&save_coloring, "-sc", "--save-coloring",
                  "-no-sc", "--dont-save-coloring",
                  "Save the mesh coloring generated when opening only a mesh.");
   args.AddOption(&portnum, "-p", "--listen-port",
                  "Specify the port number on which to accept connections.");
   args.AddOption(&secure, "-sec", "--secure-sockets",
                  "-no-sec", "--standard-sockets",
                  "Enable or disable GnuTLS secure sockets.");
   args.AddOption(&mac, "-mac", "--save-stream",
                  "-no-mac", "--dont-save-stream",
                  "In server mode, save incoming data to a file before"
                  " visualization.");
   args.AddOption(&stream_file, "-saved", "--saved-stream",
                  "Load a GLVis stream saved to a file.");
   args.AddOption(&window_w, "-ww", "--window-width",
                  "Set the window width.");
   args.AddOption(&window_h, "-wh", "--window-height",
                  "Set the window height.");
   args.AddOption(&window_title, "-wt", "--window-title",
                  "Set the window title.");
   args.AddOption(&c_plot_caption, "-c", "--plot-caption",
                  "Set the plot caption (visible when colorbar is visible).");
   args.AddOption(&font_name, "-fn", "--font",
                  "Set the font: <font-name>[-<font-size>].");
   args.AddOption(&multisample, "-ms", "--multisample",
                  "Set the multisampling mode (toggled with the 'A' key).");
   args.AddOption(&line_width, "-lw", "--line-width",
                  "Set the line width (multisampling off).");
   args.AddOption(&ms_line_width, "-mslw", "--multisample-line-width",
                  "Set the line width (multisampling on).");

   cout << endl
        << "       _/_/_/  _/      _/      _/  _/"          << endl
        << "    _/        _/      _/      _/        _/_/_/" << endl
        << "   _/  _/_/  _/      _/      _/  _/  _/_/"      << endl
        << "  _/    _/  _/        _/  _/    _/      _/_/"   << endl
        << "   _/_/_/  _/_/_/_/    _/      _/  _/_/_/"      << endl
        << endl ;

   args.Parse();
   if (!args.Good())
   {
      if (!args.Help())
      {
         args.PrintError(cout);
         cout << endl;
      }
      PrintSampleUsage(cout);
      args.PrintHelp(cout);
      return 1;
   }

   // set options
   if (mesh_file != string_none)
   {
      input |= 2;
   }
   if (sol_file != string_none)
   {
      input |= 4;
   }
   if (vec_sol_file != string_none)
   {
      sol_file = vec_sol_file;
      input |= 8;
   }
   if (gfunc_file != string_none)
   {
      sol_file = gfunc_file;
      is_gf = 255;
   }
   if (np > 0)
   {
      input |= 256;
   }
   if (arg_keys != string_none)
   {
      keys = arg_keys;
   }
   if (font_name != string_default)
   {
      SetFont(font_name);
   }
   if (multisample != GetMultisample())
   {
      SetMultisample(multisample);
   }
   if (line_width != Get_LineWidth())
   {
      Set_LineWidth(line_width);
   }
   if (ms_line_width != Get_MS_LineWidth())
   {
      Set_MS_LineWidth(ms_line_width);
   }
   if (c_plot_caption != string_none)
   {
      plot_caption = c_plot_caption;
   }

   GLVisGeometryRefiner.SetType(geom_ref_type);

   string data_type;

   // check for saved stream file
   if (stream_file != string_none)
   {
      ifstream ifs(stream_file);
      if (!ifs)
      {
         cout << "Can not open stream file: " << stream_file << endl;
         return 1;
      }
      ifs >> data_type >> ws;
      int ft = ReadStream(ifs, data_type);
      input_streams.Append(&ifs);
      StartVisualization(ft);
      return 0;
   }

   // check for script file
   if (script_file != string_none)
   {
      ifstream scr(script_file);
      if (!scr)
      {
         cout << "Can not open script: " << script_file << endl;
         return 1;
      }
      PlayScript(scr);
      return 0;
   }

   // print help for wrong input
   if (!(input == 1 || input == 3 || input == 7 || input == 11 || input == 259 ||
         (is_gf && (input == 3 || input == 259))))
   {
      cout << "Invalid combination of mesh/solution options!\n\n";
      PrintSampleUsage(cout);
      args.PrintHelp(cout);
      return 1;
   }

#ifndef MFEM_USE_GNUTLS
   if (secure)
   {
      cout << "The secure option can only be used when MFEM is compiled with"
           " GnuTLS support." << endl;
      return 1;
   }
#endif

   int childPID, viscount = 0, nproc = 1, proc = 0;

   // server mode, read the mesh and the solution from a socket
   if (input == 1)
   {
      // get rid of zombies
      if (multi_session)
      {
         signal(SIGCHLD, SIG_IGN);
      }

#ifdef MFEM_USE_GNUTLS
      GnuTLS_global_state *state = NULL;
      GnuTLS_session_params *params = NULL;
      if (secure)
      {
         state = new GnuTLS_global_state;
         // state->set_log_level(1000);
         string home_dir(getenv("HOME"));
         string server_dir = home_dir + "/.config/glvis/server/";
#ifndef MFEM_USE_GNUTLS_X509
         string pubkey  = server_dir + "pubring.gpg";
         string privkey = server_dir + "secring.gpg";
         string trustedkeys = server_dir + "trusted-clients.gpg";
#else
         string pubkey  = server_dir + "cert.pem";
         string privkey = server_dir + "key.pem";
         string trustedkeys = server_dir + "trusted-clients.pem";
#endif
         params = new GnuTLS_session_params(
            *state, pubkey.c_str(), privkey.c_str(),
            trustedkeys.c_str(), GNUTLS_SERVER);
         if (!params->status.good())
         {
            cout << "  public key   = " << pubkey << '\n'
                 << "  private key  = " << privkey << '\n'
                 << "  trusted keys = " << trustedkeys << endl;
            cout << "Error setting GLVis server parameters.\n"
                 "Generate your GLVis keys with:"
                 " bash glvis-keygen.sh [\"Your Name\"] [\"Your Email\"]"
                 << endl;
            delete params; delete state;
            return 3;
         }
      }
#endif

      const int backlog = 128;
      socketserver server(portnum, backlog);
      if (server.good())
      {
         cout << "Waiting for data on port " << portnum << " ..." << endl;
      }
      else
      {
         cout << "Server already running on port " << portnum << ".\n" << endl;
#ifdef MFEM_USE_GNUTLS
         delete params; delete state;
#endif
         return 2;
      }

      socketstream *isock;
#ifndef MFEM_USE_GNUTLS
      isock = new socketstream;
#else
      isock = secure ? new socketstream(*params) : new socketstream(false);
#endif
      while (1)
      {
         while (server.accept(*isock) < 0)
         {
#ifdef GLVIS_DEBUG
            cout << "GLVis: server.accept(...) failed." << endl;
#endif
         }

         *isock >> data_type >> ws;

         if (mac)
         {
            viscount++;
         }

         int par_data = 0;
         if (data_type == "parallel")
         {
            par_data = 1;
            np = 0;
            do
            {
               *isock >> nproc >> proc;
#ifdef GLVIS_DEBUG
               cout << "new connection: parallel " << nproc << ' ' << proc
                    << endl;
#endif
               if (np == 0)
               {
                  if (nproc <= 0)
                  {
                     cout << "Invalid number of processors: " << nproc << endl;
                     mfem_error();
                  }
                  input_streams.SetSize(nproc);
                  input_streams = NULL;
               }
               else
               {
                  if (nproc != input_streams.Size())
                  {
                     cout << "Unexpected number of processors: " << nproc
                          << ", expected: " << input_streams.Size() << endl;
                     mfem_error();
                  }
               }
               if (0 > proc || proc >= nproc)
               {
                  cout << "Invalid processor rank: " << proc
                       << ", number of processors: " << nproc << endl;
                  mfem_error();
               }
               if (input_streams[proc])
               {
                  cout << "Second connection attempt from processor rank: "
                       << proc << endl;
                  mfem_error();
               }

               input_streams[proc] = isock;
#ifndef MFEM_USE_GNUTLS
               isock = new socketstream;
#else
               isock = secure ? new socketstream(*params) :
                       new socketstream(false);
#endif
               np++;
               if (np == nproc)
               {
                  break;
               }
               // read next available socket stream
               while (server.accept(*isock) < 0)
               {
#ifdef GLVIS_DEBUG
                  cout << "GLVis: server.accept(...) failed." << endl;
#endif
               }
               *isock >> data_type >> ws; // "parallel"
               if (data_type != "parallel")
               {
                  cout << "Expected keyword \"parallel\", got \"" << data_type
                       << '"' << endl;
                  mfem_error();
               }
            }
            while (1);
         }

         char tmp_file[50];
         if (multi_session)
         {
            if (mac)
            {
               sprintf(tmp_file,"glvis-saved.%04d",viscount);
               ofstream ofs(tmp_file);
               if (!par_data)
               {
                  ofs << data_type << '\n';
                  ofs << isock->rdbuf();
                  isock->close();
               }
               else
               {
                  ReadInputStreams();
                  CloseInputStreams(false);
                  ofs.precision(8);
                  ofs << "solution\n";
                  mesh->Print(ofs);
                  grid_f->Save(ofs);
                  delete grid_f; grid_f = NULL;
                  delete mesh; mesh = NULL;
               }
               ofs.close();
               cout << "Data saved in " << tmp_file << endl;
            }
            childPID = fork();
         }
         else
         {
            childPID = 0;
         }

         switch (childPID)
         {
            case -1:
               cout << "The process couldn't fork. Exit." << endl;
               exit(1);

            case 0:                       // This is the child process
               server.close();
               if (mac)
               {
                  // exec ourself
                  const char *args[4] = { argv[0], "-saved", tmp_file, NULL };
                  execve(args[0], (char* const*)args, environ);
                  exit(0);
               }
               else
               {
                  if (multi_session)
                  {
                     signal(SIGINT, SIG_IGN);
                  }
                  int ft;
                  if (!par_data)
                  {
                     ft = ReadStream(*isock, data_type);
                     input_streams.Append(isock);
                  }
                  else
                  {
                     delete isock;
                     ft = ReadInputStreams();
                  }
                  StartVisualization(ft);
                  CloseInputStreams(false);
                  exit(0);
               }

            default :                     // This is the parent process
               if (!par_data)
               {
                  isock->rdbuf()->socketbuf::close();
               }
               else
               {
                  CloseInputStreams(true);
               }
         }
      }
#ifdef MFEM_USE_GNUTLS
      delete params; delete state;
#endif
   }
   else  // input != 1, non-server mode
   {
      if (input & 256)
      {
         ReadParallel();
      }
      else
      {
         ReadSerial();
      }

      int window_err;
      double mesh_range = -1.0;
      if (mesh->SpaceDimension() == 2)
      {
         if ((input & 8) == 0)
         {
            VisualizationSceneSolution *vss;
            window_err = InitVis((input & 4) ? 0 : 2);
            if (!window_err)
            {
               if ((input & 4) == 0)
               {
                  Set_Palette(4);   // Set_Palette(11);
               }
               vs = vss = new VisualizationSceneSolution(*mesh, sol);
               if (is_gf)
               {
                  vss->SetGridFunction(*grid_f);
               }
               if ((input & 4) == 0)
               {
                  vs->OrthogonalProjection = 1;
                  vs->light = 0;
                  vs->Zoom(1.8);
                  if (grid_f)
                  {
                     mesh_range = grid_f->Max() + 1.0;
                  }
                  else
                  {
                     mesh_range = sol.Max() + 1.0;
                  }
               }
            }
         }
         else
         {
            window_err = InitVis(1);
            if (!window_err)
            {
               if (is_gf)
               {
                  vs = new VisualizationSceneVector(*grid_f);
               }
               else
               {
                  vs = new VisualizationSceneVector(*mesh, solu, solv);
               }
            }
         }
      }
      else // 3D
      {
         if ((input & 8) == 0 && (input & 512) == 0)
         {
            VisualizationSceneSolution3d *vss;
            window_err = InitVis((input & 4) ? 0 : 2);
            if (!window_err)
            {
               vs = vss = new VisualizationSceneSolution3d(*mesh, sol);
               if (is_gf)
               {
                  vss->SetGridFunction(grid_f);
               }
               if ((input & 4) == 0)
               {
                  if (mesh->Dimension() == 3)
                  {
                     // Set_Palette(4);
                     Set_Palette(11);
                     Set_Material_And_Light(4,3);
                  }
                  else
                  {
                     Set_Palette(4);
                  }
                  vss->ToggleDrawAxes();
                  vss->ToggleDrawMesh();
                  if (grid_f)
                  {
                     mesh_range = grid_f->Max() + 1.0;
                  }
                  else
                  {
                     mesh_range = sol.Max() + 1.0;
                  }
               }
            }
         }
         else
         {
            window_err = InitVis(1);
            if (!window_err)
            {
               if (is_gf)
               {
                  grid_f = ProjectVectorFEGridFunction(grid_f);
                  vs = new VisualizationSceneVector3d(*grid_f);
               }
               else
               {
                  vs = new VisualizationSceneVector3d(*mesh, solu, solv, solw);
               }
            }
         }
      }
      if (!window_err)
      {
         // increase the refinement factors if visualizing a GridFunction
         if (is_gf)
         {
            vs->AutoRefine();
            vs->SetShading(2, true);
         }
         if (mesh_range > 0.0)
         {
            vs->SetValueRange(-mesh_range, mesh_range);
            vs->SetAutoscale(0);
         }
         if (mesh->SpaceDimension() == 2 && (input & 12) == 0)
         {
            SetVisualizationScene(vs, 2, keys.c_str());
         }
         else
         {
            SetVisualizationScene(vs, 3, keys.c_str());
         }
         KillVisualization(); // deletes vs
         if (is_gf) { delete grid_f; }
         delete mesh;
      }
      else
      {
         cerr << "Initializing the visualization failed." << endl;
         return 1;
      }
   }

   cout << "Thank you for using GLVis." << endl;

   return 0;
}


void PrintSampleUsage(ostream &out)
{
   out <<
       "Start a GLVis server:\n"
       "   glvis\n"
       "Visualize a mesh:\n"
       "   glvis -m <mesh_file>\n"
       "Visualize mesh and solution (grid function):\n"
       "   glvis -m <mesh_file> -g <grid_function_file> [-gc <component>]\n"
       "Visualize parallel mesh and solution (grid function):\n"
       "   glvis -np <#proc> -m <mesh_prefix> [-g <grid_function_prefix>]\n\n"
       "All Options:\n";
}


void ReadSerial()
{
   // get the mesh from a file
   ifgzstream meshin(mesh_file);
   if (!meshin)
   {
      cerr << "Can not open mesh file " << mesh_file << ". Exit.\n";
      exit(1);
   }

   mesh = new Mesh(meshin, 1, 0, fix_elem_orient);

   if (is_gf || (input & 4) || (input & 8))
   {
      // get the solution from file
      bool freesolin = false;
      ifgzstream *solin = NULL;
      if (!strcmp(mesh_file,sol_file))
      {
         solin = &meshin;
      }
      else
      {
         solin = new ifgzstream(sol_file);
         freesolin = true;
         if (!(*solin))
         {
            cerr << "Can not open solution file " << sol_file << ". Exit.\n";
            exit(1);
         }
      }

      if (is_gf)
      {
         grid_f = new GridFunction(mesh, *solin);
         SetGridFunction();
      }
      else if (input & 4)
      {
         // get rid of NetGen's info line
         char buff[128];
         solin->getline(buff,128);
         sol.Load(*solin, mesh->GetNV());
      }
      else if (input & 8)
      {
         solu.Load(*solin, mesh->GetNV());
         solv.Load(*solin, mesh->GetNV());
         if (mesh->SpaceDimension() == 3)
         {
            solw.Load(*solin, mesh->GetNV());
         }
      }
      if (freesolin)
      {
         delete solin;
      }
   }
   else
   {
      SetMeshSolution(mesh, grid_f, save_coloring);
   }

   Extrude1DMeshAndSolution(&mesh, &grid_f, &sol);
}


void SetGridFunction()
{
   if (gf_component != -1)
   {
      if (gf_component < 0 || gf_component >= grid_f->VectorDim())
      {
         cerr << "Invalid component " << gf_component << '.' << endl;
         exit(1);
      }
      FiniteElementSpace *ofes = grid_f->FESpace();
      FiniteElementCollection *fec =
         FiniteElementCollection::New(ofes->FEColl()->Name());
      FiniteElementSpace *fes = new FiniteElementSpace(mesh, fec);
      GridFunction *new_gf = new GridFunction(fes);
      new_gf->MakeOwner(fec);
      for (int i = 0; i < new_gf->Size(); i++)
      {
         (*new_gf)(i) = (*grid_f)(ofes->DofToVDof(i, gf_component));
      }
      delete grid_f;
      grid_f = new_gf;
   }
   if (grid_f->VectorDim() == 1)
   {
      grid_f->GetNodalValues(sol);
      input |= 4;
   }
   else
   {
      input |= 8;
   }
}


void SetMeshSolution(Mesh *mesh, GridFunction *&grid_f, bool save_coloring)
{
   if (1) // checkerboard solution
   {
      FiniteElementCollection *cfec;
      if (mesh->Dimension() == 1)
      {
         cfec = new L2_FECollection(0, 1);
      }
      else if (mesh->Dimension() == 2)
      {
         cfec = new Const2DFECollection;
      }
      else
      {
         cfec = new Const3DFECollection;
      }
      FiniteElementSpace *cfes = new FiniteElementSpace(mesh, cfec);
      grid_f = new GridFunction(cfes);
      grid_f->MakeOwner(cfec);
      {
         Array<int> coloring;
         srandom(time(0));
         double a = double(random()) / (double(RAND_MAX) + 1.);
         int el0 = (int)floor(a * mesh->GetNE());
         cout << "Generating coloring starting with element " << el0+1
              << " / " << mesh->GetNE() << endl;
         mesh->GetElementColoring(coloring, el0);
         for (int i = 0; i < coloring.Size(); i++)
         {
            (*grid_f)(i) = coloring[i];
         }
         cout << "Number of colors: " << grid_f->Max() + 1 << endl;
      }
      grid_f->GetNodalValues(sol);
      is_gf = 1;
      if (save_coloring)
      {
         const char col_fname[] = "GLVis_coloring.gf";
         ofstream fgrid(col_fname);
         cout << "Saving the coloring function -> " << flush;
         grid_f->Save(fgrid);
         cout << col_fname << endl;
      }
   }
   else // zero solution
   {
      sol.SetSize (mesh -> GetNV());
      sol = 0.0;
   }
}


void ReadParallel()
{
   int err;

   if (is_gf)
   {
      err = ReadParMeshAndGridFunction(np, mesh_file, sol_file,
                                       &mesh, &grid_f, keep_attr);
      if (!err)
      {
         SetGridFunction();
      }
   }
   else
   {
      err = ReadParMeshAndGridFunction(np, mesh_file, NULL,
                                       &mesh, NULL, keep_attr);
      if (!err)
      {
         SetMeshSolution(mesh, grid_f, save_coloring);
      }
   }

   if (err)
   {
      exit(1);
   }

   Extrude1DMeshAndSolution(&mesh, &grid_f, &sol);
}

int ReadParMeshAndGridFunction(int np, const char *mesh_prefix,
                               const char *sol_prefix, Mesh **mesh_p,
                               GridFunction **sol_p, int keep_attr)
{
   Array<Mesh *> mesh_array;
   Array<named_ifgzstream *> meshfiles;

   mesh_array.SetSize(np);
   meshfiles.SetSize(np);
   for (int p = 0; p < np; p++)
   {
      ostringstream fname;
      fname << mesh_prefix << '.' << setfill('0') << setw(pad_digits) << p;
      meshfiles[p] = new named_ifgzstream(fname.str().c_str());
      if (!(*meshfiles[p]))
      {
         cerr << "Can not open mesh file: " << fname.str().c_str()
              << '!' << endl;
         for (p--; p >= 0; p--)
         {
            delete mesh_array[p];
            delete meshfiles[p];
         }
         return 1;
      }
      mesh_array[p] = new Mesh(*meshfiles[p], 1, 0, fix_elem_orient);
      if (!keep_attr)
      {
         // set element and boundary attributes to be the processor number + 1
         for (int i = 0; i < mesh_array[p]->GetNE(); i++)
         {
            mesh_array[p]->GetElement(i)->SetAttribute(p+1);
         }
         for (int i = 0; i < mesh_array[p]->GetNBE(); i++)
         {
            mesh_array[p]->GetBdrElement(i)->SetAttribute(p+1);
         }
      }
   }
   *mesh_p = new Mesh(mesh_array, np);

   if (sol_prefix && sol_p)
   {
      Array<GridFunction *> gf_array(np);
      bool samefile = !strcmp(sol_prefix,mesh_prefix);
      for (int p = 0; p < np; p++)
      {
         if (!samefile)
         {
            ostringstream fname;
            fname << sol_prefix << '.' << setfill('0') << setw(pad_digits) << p;
            ifgzstream solfile(fname.str().c_str());
            if (!solfile)
            {
               cerr << "Can not open solution file " << fname.str().c_str()
                    << '!' << endl;
               for (p--; p >= 0; p--)
               {
                  delete gf_array[p];
               }
               delete *mesh_p;
               *mesh_p = NULL;
               for (p = 0; p < np; p++)
               {
                  delete mesh_array[np-1-p];
                  delete meshfiles[np-1-p];
               }
               return 2;
            }
            gf_array[p] = new GridFunction(mesh_array[p], solfile);
         }
         else  // mesh and solution in the same file
         {
            gf_array[p] = new GridFunction(mesh_array[p], *meshfiles[p]);
         }
      }
      *sol_p = new GridFunction(*mesh_p, gf_array, np);

      for (int p = 0; p < np; p++)
      {
         delete gf_array[np-1-p];
      }
   }

   for (int p = 0; p < np; p++)
   {
      delete mesh_array[np-1-p];
      delete meshfiles[np-1-p];
   }

   return 0;
}

int ReadInputStreams()
{
   int nproc = input_streams.Size();
   Array<Mesh *> mesh_array(nproc);
   Array<GridFunction *> gf_array(nproc);
   string data_type;

   int gf_count = 0;
   int field_type = 0;

   for (int p = 0; p < nproc; p++)
   {
#ifdef GLVIS_DEBUG
      cout << "connection[" << p << "]: reading initial data ... " << flush;
#endif
      istream &isock = *input_streams[p];
      // assuming the "parallel nproc p" part of the stream has been read
      isock >> ws >> data_type >> ws; // "*_data" / "mesh" / "solution"
#ifdef GLVIS_DEBUG
      cout << " type " << data_type << " ... " << flush;
#endif
      mesh_array[p] = new Mesh(isock, 1, 0, fix_elem_orient);
      if (!keep_attr)
      {
         // set element and boundary attributes to proc+1
         for (int i = 0; i < mesh_array[p]->GetNE(); i++)
         {
            mesh_array[p]->GetElement(i)->SetAttribute(p+1);
         }
         for (int i = 0; i < mesh_array[p]->GetNBE(); i++)
         {
            mesh_array[p]->GetBdrElement(i)->SetAttribute(p+1);
         }
      }
      gf_array[p] = NULL;
      if (data_type != "mesh")
      {
         gf_array[p] = new GridFunction(mesh_array[p], isock);
         gf_count++;
      }
#ifdef GLVIS_DEBUG
      cout << "done." << endl;
#endif
   }

   if (gf_count > 0 && gf_count != nproc)
   {
      mfem_error("Input streams contain a mixture of data types!");
   }

   mesh = new Mesh(mesh_array, nproc);
   if (gf_count == 0)
   {
      SetMeshSolution(mesh, grid_f, save_coloring);
      field_type = 2;
   }
   else
   {
      grid_f = new GridFunction(mesh, gf_array, nproc);
      field_type = (grid_f->VectorDim() == 1) ? 0 : 1;
   }

   for (int p = 0; p < nproc; p++)
   {
      delete mesh_array[nproc-1-p];
      delete gf_array[nproc-1-p];
   }

   Extrude1DMeshAndSolution(&mesh, &grid_f, NULL);

   return field_type;
}

void CloseInputStreams(bool parent)
{
   for (int i = 0; i < input_streams.Size(); i++)
   {
      if (parent)
      {
         socketstream *sock = dynamic_cast<socketstream*>(input_streams[i]);
         if (sock) { sock->rdbuf()->socketbuf::close(); }
      }
      delete input_streams[i];
   }
   input_streams.DeleteAll();
}

// 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.
GridFunction *ProjectVectorFEGridFunction(GridFunction *gf)
{
   if ((gf->VectorDim() == 3) && (gf->FESpace()->GetVDim() == 1))
   {
      int p = gf->FESpace()->GetOrder(0);
      cout << "Switching to order " << p
           << " discontinuous vector grid function..." << endl;
      int dim = gf->FESpace()->GetMesh()->Dimension();
      FiniteElementCollection *d_fec = new L2_FECollection(p, dim, 1);
      FiniteElementSpace *d_fespace =
         new FiniteElementSpace(gf->FESpace()->GetMesh(), d_fec, 3);
      GridFunction *d_gf = new GridFunction(d_fespace);
      d_gf->MakeOwner(d_fec);
      gf->ProjectVectorFieldOn(*d_gf);
      delete gf;
      return d_gf;
   }
   return gf;
}

void Extrude1DMeshAndSolution(Mesh **mesh_p, GridFunction **grid_f_p,
                              Vector *sol)
{
   Mesh *mesh = *mesh_p;

   if (mesh->Dimension() != 1 || mesh->SpaceDimension() != 1)
   {
      return;
   }

   // find xmin and xmax over the vertices of the 1D mesh
   double xmin = numeric_limits<double>::infinity();
   double xmax = -xmin;
   for (int i = 0; i < mesh->GetNV(); i++)
   {
      const double x = mesh->GetVertex(i)[0];
      if (x < xmin)
      {
         xmin = x;
      }
      if (x > xmax)
      {
         xmax = x;
      }
   }

   Mesh *mesh2d = Extrude1D(mesh, 1, 0.1*(xmax - xmin));

   if (grid_f_p && *grid_f_p)
   {
      GridFunction *grid_f_2d =
         Extrude1DGridFunction(mesh, mesh2d, *grid_f_p, 1);

      delete *grid_f_p;
      *grid_f_p = grid_f_2d;
   }
   if (sol && sol->Size() == mesh->GetNV())
   {
      Vector sol2d(mesh2d->GetNV());
      for (int i = 0; i < mesh->GetNV(); i++)
      {
         sol2d(2*i+0) = sol2d(2*i+1) = (*sol)(i);
      }
      *sol = sol2d;
   }

   delete mesh;
   *mesh_p = mesh2d;
}