From 832942340400fbc20a79c095471fbfb2c79ab993 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 11:28:51 -0800 Subject: [PATCH 01/44] Added visualization of complex functions from file. --- lib/data_state.cpp | 121 ++++++++++++++++++++++++++++++++++++++++++++ lib/data_state.hpp | 37 ++++++++++++++ lib/file_reader.cpp | 22 +++++++- lib/window.cpp | 20 +++++++- lib/window.hpp | 6 +++ 5 files changed, 203 insertions(+), 3 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index 461f233b..914dd384 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -40,6 +40,7 @@ DataState &DataState::operator=(DataState &&ss) internal = std::move(ss.internal); type = ss.type; + cmplx_sol = ss.cmplx_sol; quad_sol = ss.quad_sol; sol = std::move(ss.sol); @@ -67,6 +68,7 @@ void DataState::SetMesh(std::unique_ptr &&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); } } @@ -80,11 +82,41 @@ void DataState::SetGridFunction( std::unique_ptr &&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); } +void DataState::SetCmplxGridFunction(mfem::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 + &&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::SetQuadFunction(mfem::QuadratureFunction *qf, int component) { if (quad_f.get() != qf) @@ -93,6 +125,7 @@ void DataState::SetQuadFunction(mfem::QuadratureFunction *qf, int component) quad_sol = QuadSolution::NONE; } internal.quad_f.reset(qf); + internal.cgrid_f.reset(); SetQuadFunctionSolution(component); } @@ -105,6 +138,7 @@ void DataState::SetQuadFunction( quad_sol = QuadSolution::NONE; } internal.quad_f = std::move(pqf); + internal.cgrid_f.reset(); SetQuadFunctionSolution(component); } @@ -369,6 +403,48 @@ void DataState::SetGridFunctionSolution(int gf_component) } } +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) @@ -409,6 +485,51 @@ void DataState::SetQuadFunctionSolution(int qf_component) SetQuadSolution(); } +void DataState::SetComplexSolution(ComplexSolution cmplx_type) +{ + double (*cmplx_2_scalar)(double, double); + const char *str_cmplx_2_scalar; + + switch (cmplx_type) + { + case ComplexSolution::Magnitude: + cmplx_2_scalar = [](double r, double i) { return hypot(r,i); }; + str_cmplx_2_scalar = "magnitude"; + break; + case ComplexSolution::Phase: + cmplx_2_scalar = [](double r, double i) { return atan2(r,i); }; + str_cmplx_2_scalar = "phase"; + break; + case ComplexSolution::Real: + cmplx_2_scalar = [](double r, double i) { return r; }; + str_cmplx_2_scalar = "real part"; + break; + case ComplexSolution::Imag: + cmplx_2_scalar = [](double r, double i) { return i; }; + str_cmplx_2_scalar = "imaginary part"; + break; + default: + cout << "Unknown complex data representation" << endl; + return; + } + + cout << "Representing complex function by: " << str_cmplx_2_scalar << endl; + GridFunction *gf = new GridFunction(cgrid_f->FESpace()); + internal.grid_f.reset(gf); + for (int i = 0; i < gf->Size(); i++) + { + (*gf)(i) = cmplx_2_scalar(cgrid_f->real()(i), cgrid_f->imag()(i)); + } + + cmplx_sol = cmplx_type; +} + +void DataState::SwitchComplexSolution(ComplexSolution cmplx_type) +{ + SetComplexSolution(cmplx_type); + ExtrudeMeshAndSolution(); +} + void DataState::SetQuadSolution(QuadSolution quad_type) { // assume identical order diff --git a/lib/data_state.hpp b/lib/data_state.hpp index 3ee65235..031dfbd2 100644 --- a/lib/data_state.hpp +++ b/lib/data_state.hpp @@ -44,6 +44,19 @@ struct DataState MAX }; + enum class ComplexSolution + { + NONE = -1, + MIN = -1, + //---------- + Magnitude, + Phase, + Real, + Imag, + //---------- + MAX + }; + private: friend class StreamReader; friend class FileReader; @@ -52,14 +65,17 @@ struct DataState std::unique_ptr mesh; std::unique_ptr mesh_quad; std::unique_ptr grid_f; + std::unique_ptr cgrid_f; std::unique_ptr quad_f; std::unique_ptr data_coll; } 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); public: @@ -67,6 +83,7 @@ struct DataState const std::unique_ptr &mesh{internal.mesh}; const std::unique_ptr &mesh_quad{internal.mesh_quad}; const std::unique_ptr &grid_f{internal.grid_f}; + const std::unique_ptr &cgrid_f{internal.cgrid_f}; const std::unique_ptr &quad_f{internal.quad_f}; const std::unique_ptr &data_coll{internal.data_coll}; @@ -102,6 +119,17 @@ struct DataState void SetGridFunction(std::unique_ptr &&pgf, int component = -1); + /// Set a complex grid function (plain pointer version) + /** Note that ownership is passed from the caller. + @see SetCmplxGridFunction(std::unique_ptr &&, 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 &&pgf, + int component = -1); + /// Set a quadrature function (plain pointer version) /** Note that ownership is passed from the caller. @see SetQuadFunction(std::unique_ptr &&, int ) */ @@ -141,6 +169,15 @@ 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); + + /// Switch the complex function representation + void SwitchComplexSolution(ComplexSolution type); + + /// 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); diff --git a/lib/file_reader.cpp b/lib/file_reader.cpp index 7efc897e..6346abe0 100644 --- a/lib/file_reader.cpp +++ b/lib/file_reader.cpp @@ -12,6 +12,7 @@ #include "file_reader.hpp" #include +#include using namespace std; using namespace mfem; @@ -54,8 +55,25 @@ 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; + { + string buff; + auto pos = solin->tellg(); + *solin >> ws; + getline(*solin, buff); + mfem::filter_dos(buff); + const bool cmplx = (buff == "ComplexGridFunction"); + solin->seekg(pos); + if (cmplx) + { + 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); diff --git a/lib/window.cpp b/lib/window.cpp index 78ad78a0..e7c7ce6e 100644 --- a/lib/window.cpp +++ b/lib/window.cpp @@ -73,7 +73,11 @@ bool Window::GLVisInitVis(StreamCollection input_streams) locwin = this; - if (data_state.quad_f) + if (data_state.cgrid_f) + { + wnd->setOnKeyDown('Q', SwitchComplexSolution); + } + else if (data_state.quad_f) { wnd->setOnKeyDown('Q', SwitchQuadSolution); } @@ -194,6 +198,12 @@ void Window::GLVisStartVis() std::cout << "GLVis window closed." << std::endl; } +void Window::SwitchComplexSolution(DataState::ComplexSolution cmplx_type) +{ + data_state.SwitchComplexSolution(cmplx_type); + ResetMeshAndSolution(data_state); +} + void Window::SwitchQuadSolution(DataState::QuadSolution quad_type) { data_state.SwitchQuadSolution(quad_type); @@ -260,6 +270,14 @@ void Window::ResetMeshAndSolution(DataState &ss) thread_local Window *Window::locwin = NULL; +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 5bea004b..b742e5aa 100644 --- a/lib/window.hpp +++ b/lib/window.hpp @@ -61,6 +61,9 @@ 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); @@ -82,6 +85,9 @@ struct Window /// 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(); }; From 47574873401c461f3c8a6b431a4c369241e52168 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 13:15:07 -0800 Subject: [PATCH 02/44] Added support for vector spaces. --- lib/data_state.cpp | 87 ++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 80 insertions(+), 7 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index 914dd384..ec488485 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -487,8 +487,9 @@ void DataState::SetQuadFunctionSolution(int qf_component) void DataState::SetComplexSolution(ComplexSolution cmplx_type) { - double (*cmplx_2_scalar)(double, double); + double (*cmplx_2_scalar)(double, double) = NULL; const char *str_cmplx_2_scalar; + GridFunction *gf = NULL; switch (cmplx_type) { @@ -501,11 +502,11 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type) str_cmplx_2_scalar = "phase"; break; case ComplexSolution::Real: - cmplx_2_scalar = [](double r, double i) { return r; }; + gf = new GridFunction(cgrid_f->FESpace(), cgrid_f->real(), 0); str_cmplx_2_scalar = "real part"; break; case ComplexSolution::Imag: - cmplx_2_scalar = [](double r, double i) { return i; }; + gf = new GridFunction(cgrid_f->FESpace(), cgrid_f->imag(), 0); str_cmplx_2_scalar = "imaginary part"; break; default: @@ -514,12 +515,84 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type) } cout << "Representing complex function by: " << str_cmplx_2_scalar << endl; - GridFunction *gf = new GridFunction(cgrid_f->FESpace()); - internal.grid_f.reset(gf); - for (int i = 0; i < gf->Size(); i++) + if (!gf) { - (*gf)(i) = cmplx_2_scalar(cgrid_f->real()(i), cgrid_f->imag()(i)); + gf = new GridFunction(cgrid_f->FESpace()); + const FiniteElementSpace *fes = cgrid_f->FESpace(); + const Mesh *mesh = fes->GetMesh(); + const int dim = mesh->Dimension(); + if (cgrid_f->FESpace()->FEColl()->GetMapType(dim) == FiniteElement::VALUE) + { + for (int i = 0; i < gf->Size(); i++) + { + (*gf)(i) = cmplx_2_scalar(cgrid_f->real()(i), cgrid_f->imag()(i)); + } + } + else + { + const int vdim = fes->GetVDim(); + const int sdim = fes->GetVectorDim(); + Array vdofs; + ElementTransformation *Tr; + Vector r_data, i_data, z_data; + DenseMatrix r_vals, i_vals, z_vals; + Vector r_vec, i_vec, z_vec; + 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); + z_data.SetSize(vdofs.Size()); + + Tr = fes->GetElementTransformation(z); + const FiniteElement *fe = fes->GetFE(z); + const IntegrationRule &ir = fe->GetNodes(); + const int nnp = ir.GetNPoints(); + r_vals.Reset(r_data.GetData(), nnp, vdim); + i_vals.Reset(i_data.GetData(), nnp, vdim); + z_vals.Reset(z_data.GetData(), nnp, vdim); + + if (fe->GetRangeType() == FiniteElement::SCALAR) + { + shape.SetSize(nnp); + } + else + { + vshape.SetSize(nnp, sdim); + } + for (int n = 0; n < nnp; n++) + { + const IntegrationPoint &ip = ir.IntPoint(n); + Tr->SetIntPoint(&ip); + double w; + if (fe->GetRangeType() == FiniteElement::SCALAR) + { + fe->CalcPhysShape(*Tr, shape); + w = shape(n); + } + else + { + fe->CalcPhysVShape(*Tr, vshape); + Vector vec(sdim); + vshape.GetRow(n, vec); + w = vec.Norml2(); + } + for (int d = 0; d < vdim; d++) + { + const double rval = r_vals(n,d) * w; + const double ival = i_vals(n,d) * w; + const double zval = cmplx_2_scalar(rval, ival); + z_vals(n,d) = (w != 0.)?(zval / w):(0.); + } + } + + gf->SetSubVector(vdofs, z_data); + } + } } + internal.grid_f.reset(gf); cmplx_sol = cmplx_type; } From c306818999e6de707422af8fed07097730c532f8 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 14:10:48 -0800 Subject: [PATCH 03/44] Added loading of parallel complex grid function. --- lib/file_reader.cpp | 83 ++++++++++++++++++++++++++++++++++++++------- lib/file_reader.hpp | 2 ++ 2 files changed, 73 insertions(+), 12 deletions(-) diff --git a/lib/file_reader.cpp b/lib/file_reader.cpp index 6346abe0..6a81c149 100644 --- a/lib/file_reader.cpp +++ b/lib/file_reader.cpp @@ -56,14 +56,7 @@ int FileReader::ReadSerial(FileReader::FileType ft, const char *mesh_file, { case FileType::GRID_FUNC: { - string buff; - auto pos = solin->tellg(); - *solin >> ws; - getline(*solin, buff); - mfem::filter_dos(buff); - const bool cmplx = (buff == "ComplexGridFunction"); - solin->seekg(pos); - if (cmplx) + if (CheckStreamIsComplex(*solin)) { data.SetCmplxGridFunction(new ComplexGridFunction(data.mesh.get(), *solin), component); @@ -135,6 +128,19 @@ int FileReader::ReadParallel(int np, FileType ft, const char *mesh_file, return read_err; } +bool FileReader::CheckStreamIsComplex(std::istream &solin, bool parallel) const +{ + 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) { @@ -150,6 +156,7 @@ int FileReader::ReadParMeshAndGridFunction(int np, const char *mesh_prefix, std::vector mesh_array(np); std::vector gf_array(np); + std::vector cgf_array(np); int read_err = 0; for (int p = 0; p < np; p++) @@ -195,23 +202,74 @@ 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); + } + else + { + 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], meshfile); + if (CheckStreamIsComplex(meshfile)) + { + cgf_array[p] = new ComplexGridFunction(mesh_array[p], meshfile); + } + else + { + gf_array[p] = new GridFunction(mesh_array[p], meshfile); + } } } } + if (!read_err) + { + bool bcmplx{}, breal{}; + for (int p = 0; p < np; p++) + { + if (cgf_array[p]) { bcmplx = true; } + if (gf_array[p]) { breal = true; } + } + if (bcmplx && breal) + { + cerr << "Inconsistent input files" << 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(new GridFunction(data.mesh.get(), gf_array.data(), np), - component); + if (cgf_array[0]) + { + std::vector 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(data.mesh.get(), r_array.data(), np); + GridFunction *igf = new GridFunction(data.mesh.get(), i_array.data(), np); + ComplexGridFunction *cgf = new ComplexGridFunction(rgf->FESpace()); + cgf->MakeOwner(rgf->OwnFEC()); + rgf->MakeOwner(NULL); + cgf->real() = *rgf; + cgf->imag() = *igf; + delete rgf; + delete igf; + data.SetCmplxGridFunction(cgf, component); + } + else + data.SetGridFunction(new GridFunction(data.mesh.get(), gf_array.data(), np), + component); } else { @@ -222,6 +280,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 762cd1bc..337704cd 100644 --- a/lib/file_reader.hpp +++ b/lib/file_reader.hpp @@ -21,6 +21,8 @@ class FileReader DataState &data; int pad_digits; + bool CheckStreamIsComplex(std::istream &sol, bool parallel = false) const; + int ReadParMeshAndGridFunction(int np, const char *mesh_prefix, const char *sol_prefix, int component = -1); int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, From d0eee5784f4471db3f44ea978e449246026d52f7 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 14:26:45 -0800 Subject: [PATCH 04/44] Added streaming of complex solution. --- lib/stream_reader.cpp | 6 ++++++ lib/threads.cpp | 27 ++++++++++++++++++++------- 2 files changed, 26 insertions(+), 7 deletions(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 24b40d85..2fd5d94c 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -21,6 +21,7 @@ enum class Command { Mesh, Solution, + ComplexSolution, Quadrature, Fem2D, VFem2D, @@ -68,6 +69,7 @@ StreamCommands::StreamCommands() { (*this)[Command::Mesh] = {"mesh", false, "", "Visualize the mesh."}; (*this)[Command::Solution] = {"solution", false, " ", "Visualize the solution."}; + (*this)[Command::ComplexSolution] = {"csolution", false, " ", "Visualize the complex solution."}; (*this)[Command::Quadrature] = {"quadrature", false, " ", "Visualize the quadrature."}; (*this)[Command::Fem2D] = {"fem2d_data", false, " ", "Visualize the 2D scalar data."}; (*this)[Command::VFem2D] = {"vfem2d_data", false, " ", "Visualize the 2D vector data."}; @@ -156,6 +158,10 @@ int StreamReader::ReadStream( data.SetMesh(new Mesh(is, 1, 0, data.fix_elem_orient)); data.SetGridFunction(new GridFunction(data.mesh.get(), is)); break; + case Command::ComplexSolution: + data.SetMesh(new Mesh(is, 1, 0, data.fix_elem_orient)); + data.SetCmplxGridFunction(new ComplexGridFunction(data.mesh.get(), is)); + break; case Command::Quadrature: data.SetMesh(new Mesh(is, 1, 0, data.fix_elem_orient)); data.SetQuadFunction(new QuadratureFunction(data.mesh.get(), is)); diff --git a/lib/threads.cpp b/lib/threads.cpp index c53717cf..66634df5 100644 --- a/lib/threads.cpp +++ b/lib/threads.cpp @@ -469,13 +469,7 @@ int GLVisCommand::Execute() double mesh_range = -1.0; if (!new_state.grid_f) { - if (!new_state.quad_f) - { - new_state.save_coloring = false; - new_state.SetMeshSolution(); - mesh_range = new_state.grid_f->Max() + 1.0; - } - else + if (new_state.quad_f) { auto qs = win.data_state.GetQuadSolution(); if (qs != DataState::QuadSolution::NONE) @@ -488,6 +482,25 @@ 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(); + } + else + { + new_state.save_coloring = false; + new_state.SetMeshSolution(); + mesh_range = new_state.grid_f->Max() + 1.0; + } } if (win.SetNewMeshAndSolution(std::move(new_state))) { From 437aa1c6cd41137533a665613d4b7912c1e287cf Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 16:22:48 -0800 Subject: [PATCH 05/44] Fixed setting of the representation state during updates. --- lib/threads.cpp | 49 +++++++++++++++++++++++-------------------------- 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/lib/threads.cpp b/lib/threads.cpp index 66634df5..a0aa6509 100644 --- a/lib/threads.cpp +++ b/lib/threads.cpp @@ -467,40 +467,37 @@ int GLVisCommand::Execute() case NEW_MESH_AND_SOLUTION: { double mesh_range = -1.0; - if (!new_state.grid_f) + if (new_state.quad_f) { - if (new_state.quad_f) + auto qs = win.data_state.GetQuadSolution(); + if (qs != DataState::QuadSolution::NONE) { - auto qs = win.data_state.GetQuadSolution(); - if (qs != DataState::QuadSolution::NONE) - { - new_state.SetQuadSolution(qs); - } - else - { - new_state.SetQuadSolution(); - } - new_state.ExtrudeMeshAndSolution(); + new_state.SetQuadSolution(qs); } - else if (new_state.cgrid_f) + else { - auto cs = win.data_state.GetComplexSolution(); - if (cs != DataState::ComplexSolution::NONE) - { - new_state.SetComplexSolution(cs); - } - else - { - new_state.SetComplexSolution(); - } - new_state.ExtrudeMeshAndSolution(); + new_state.SetQuadSolution(); + } + 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.save_coloring = false; - new_state.SetMeshSolution(); - mesh_range = new_state.grid_f->Max() + 1.0; + new_state.SetComplexSolution(); } + new_state.ExtrudeMeshAndSolution(); + } + else if (!new_state.grid_f) + { + new_state.save_coloring = false; + new_state.SetMeshSolution(); + mesh_range = new_state.grid_f->Max() + 1.0; } if (win.SetNewMeshAndSolution(std::move(new_state))) { From 4dbc1965e7ce0ef52ac226978d44831c3455060f Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 16:27:57 -0800 Subject: [PATCH 06/44] Added stream writer for complex gfs. --- lib/stream_reader.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 2fd5d94c..d4fafa27 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -398,6 +398,12 @@ void StreamReader::WriteStream(std::ostream &os) } data.quad_f->Save(os); } + else if (data.cgrid_f) + { + os << "csolution\n"; + data.mesh->Print(os); + data.cgrid_f->Save(os); + } else if (data.grid_f) { os << "solution\n"; From f554774ee9724379108e2de5ef97ad50b420d380 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 16:54:55 -0800 Subject: [PATCH 07/44] Added loading of parallel streams of complex gfs. --- lib/data_state.cpp | 27 +++++++++++++++++++++++++-- lib/data_state.hpp | 6 ++++++ lib/file_reader.cpp | 17 +---------------- lib/stream_reader.cpp | 16 +++++++++++++++- 4 files changed, 47 insertions(+), 19 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index ec488485..b69ab16f 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -103,8 +103,8 @@ void DataState::SetCmplxGridFunction(mfem::ComplexGridFunction *gf, SetComplexFunctionSolution(component); } -void DataState::SetCmplxGridFunction(std::unique_ptr - &&pgf, int component) +void DataState::SetCmplxGridFunction( + std::unique_ptr &&pgf, int component) { if (cgrid_f.get() != pgf.get()) { @@ -117,6 +117,29 @@ void DataState::SetCmplxGridFunction(std::unique_ptr SetComplexFunctionSolution(component); } +void DataState::SetCmplxGridFunction( + const std::vector &cgf_array, int component) +{ + const int np = cgf_array.size(); + std::vector 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); +} + void DataState::SetQuadFunction(mfem::QuadratureFunction *qf, int component) { if (quad_f.get() != qf) diff --git a/lib/data_state.hpp b/lib/data_state.hpp index 031dfbd2..20ac6efa 100644 --- a/lib/data_state.hpp +++ b/lib/data_state.hpp @@ -130,6 +130,12 @@ struct DataState void SetCmplxGridFunction(std::unique_ptr &&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 + &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 &&, int ) */ diff --git a/lib/file_reader.cpp b/lib/file_reader.cpp index 6a81c149..2197e3f0 100644 --- a/lib/file_reader.cpp +++ b/lib/file_reader.cpp @@ -250,22 +250,7 @@ int FileReader::ReadParMeshAndGridFunction(int np, const char *mesh_prefix, { if (cgf_array[0]) { - std::vector 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(data.mesh.get(), r_array.data(), np); - GridFunction *igf = new GridFunction(data.mesh.get(), i_array.data(), np); - ComplexGridFunction *cgf = new ComplexGridFunction(rgf->FESpace()); - cgf->MakeOwner(rgf->OwnFEC()); - rgf->MakeOwner(NULL); - cgf->real() = *rgf; - cgf->imag() = *igf; - delete rgf; - delete igf; - data.SetCmplxGridFunction(cgf, component); + data.SetCmplxGridFunction(cgf_array, component); } else data.SetGridFunction(new GridFunction(data.mesh.get(), gf_array.data(), np), diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index d4fafa27..5d5196bb 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -304,11 +304,13 @@ int StreamReader::ReadStreams(const StreamCollection &input_streams) const int nproc = input_streams.size(); std::vector mesh_array(nproc); std::vector gf_array(nproc); + std::vector cgf_array(nproc); std::vector 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++) @@ -336,7 +338,14 @@ int StreamReader::ReadStreams(const StreamCollection &input_streams) } } gf_array[p] = NULL; - if (data_type == "quadrature") + if (data_type == "csolution") + { + isock >> ws; + isock.ignore(3);// ignore 'Par' prefix to load as serial + cgf_array[p] = new ComplexGridFunction(mesh_array[p], isock); + cgf_count++; + } + else if (data_type == "quadrature") { qf_array[p] = new QuadratureFunction(mesh_array[p], isock); qf_count++; @@ -352,6 +361,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!"); @@ -362,6 +372,10 @@ int StreamReader::ReadStreams(const StreamCollection &input_streams) { data.SetGridFunction(new GridFunction(data.mesh.get(), gf_array.data(), nproc)); } + else if (cgf_count > 0) + { + data.SetCmplxGridFunction(cgf_array); + } else if (qf_count > 0) { data.SetQuadFunction(qf_array); From cf8d775a8c75bc4f8e71146bef490558552e7985 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 16:56:01 -0800 Subject: [PATCH 08/44] Made the file reader consistency check more strict. --- lib/file_reader.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/lib/file_reader.cpp b/lib/file_reader.cpp index 2197e3f0..27d5a23c 100644 --- a/lib/file_reader.cpp +++ b/lib/file_reader.cpp @@ -158,6 +158,9 @@ int FileReader::ReadParMeshAndGridFunction(int np, const char *mesh_prefix, std::vector gf_array(np); std::vector cgf_array(np); + int gf_count = 0; + int cgf_count = 0; + int read_err = 0; for (int p = 0; p < np; p++) { @@ -207,10 +210,12 @@ int FileReader::ReadParMeshAndGridFunction(int np, const char *mesh_prefix, 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 @@ -218,10 +223,12 @@ int FileReader::ReadParMeshAndGridFunction(int np, const char *mesh_prefix, 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++; } } } @@ -229,15 +236,10 @@ int FileReader::ReadParMeshAndGridFunction(int np, const char *mesh_prefix, if (!read_err) { - bool bcmplx{}, breal{}; - for (int p = 0; p < np; p++) - { - if (cgf_array[p]) { bcmplx = true; } - if (gf_array[p]) { breal = true; } - } - if (bcmplx && breal) + if ((gf_count > 0 && gf_count != np) + || (cgf_count > 0 && cgf_count != np)) { - cerr << "Inconsistent input files" << endl; + cerr << "Input files contain a mixture of data types!" << endl; read_err = 3; } } From 7119a69f1f75045ebb95e565bd26694e5f54b028 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 18:42:51 -0800 Subject: [PATCH 09/44] Fixed shadowing in DataState. --- lib/data_state.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index b69ab16f..db999028 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -542,8 +542,8 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type) { gf = new GridFunction(cgrid_f->FESpace()); const FiniteElementSpace *fes = cgrid_f->FESpace(); - const Mesh *mesh = fes->GetMesh(); - const int dim = mesh->Dimension(); + const Mesh *msh = fes->GetMesh(); + const int dim = msh->Dimension(); if (cgrid_f->FESpace()->FEColl()->GetMapType(dim) == FiniteElement::VALUE) { for (int i = 0; i < gf->Size(); i++) @@ -562,7 +562,7 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type) Vector r_vec, i_vec, z_vec; DenseMatrix vshape; Vector shape; - for (int z = 0; z < mesh->GetNE(); z++) + for (int z = 0; z < msh->GetNE(); z++) { fes->GetElementVDofs(z, vdofs); cgrid_f->real().GetSubVector(vdofs, r_data); From 19f336e26314c5ff9e07b8f06cd155722a271099 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 21:56:31 -0800 Subject: [PATCH 10/44] Added animation of complex phase. --- lib/aux_vis.cpp | 55 +++++++++++++++- lib/aux_vis.hpp | 3 + lib/data_state.cpp | 152 +++++++++++++++++++++++++-------------------- lib/data_state.hpp | 6 +- lib/openglvis.hpp | 2 + lib/vsdata.hpp | 2 + lib/window.cpp | 7 +++ lib/window.hpp | 3 + 8 files changed, 157 insertions(+), 73 deletions(-) diff --git a/lib/aux_vis.cpp b/lib/aux_vis.cpp index 4a260e03..3380441a 100644 --- a/lib/aux_vis.cpp +++ b/lib/aux_vis.cpp @@ -229,6 +229,9 @@ GLWindow* InitVisualization(const char name[], int x, int y, int w, int h, wnd->setOnKeyDown (SDLK_PERIOD, KeyDeletePressed); wnd->setOnKeyDown (SDLK_RETURN, KeyEnterPressed); + wnd->setOnKeyDown (SDLK_SEMICOLON, KeySemicolonPressed); + wnd->setOnKeyDown (SDLK_COLON, KeyColonPressed); + wnd->setOnKeyDown (SDLK_0, Key0Pressed); wnd->setOnKeyDown (SDLK_1, Key1Pressed); wnd->setOnKeyDown (SDLK_2, Key2Pressed); @@ -612,6 +615,9 @@ 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.01; void MainLoop() { @@ -621,13 +627,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) @@ -1436,6 +1448,45 @@ void KeyEnterPressed() CheckSpin(); } +void CheckPhaseAnim() +{ + if (fabs(phase_rate) < 1.e-2) + { + phase_rate = 0.; + } + if (phase_rate != 0.) + { + phase_anim = 1; + AddIdleFunc(MainLoop); + } + else + { + phase_anim = 0; + RemoveIdleFunc(MainLoop); + } + cout << "Phase rate: " << phase_rate << " radians / frame" << endl; +} + +void KeySemicolonPressed() +{ + if (!phase_anim) + { + phase_rate = 0.; + } + phase_rate += phase_step; + CheckPhaseAnim(); +} + +void KeyColonPressed() +{ + if (!phase_anim) + { + phase_rate = 0.; + } + phase_rate -= phase_step; + CheckPhaseAnim(); +} + void Key7Pressed() { locscene->PreRotate(1.0, 0.0, -1.0, 0.0); diff --git a/lib/aux_vis.hpp b/lib/aux_vis.hpp index f0e1916d..c77fc2b7 100644 --- a/lib/aux_vis.hpp +++ b/lib/aux_vis.hpp @@ -88,6 +88,9 @@ void Key0Pressed(); void KeyDeletePressed(); void KeyEnterPressed(); +void KeySemicolonPressed(); +void KeyColonPressed(); + void KeyLeftPressed(GLenum); void KeyRightPressed(GLenum); void KeyUpPressed(GLenum); diff --git a/lib/data_state.cpp b/lib/data_state.cpp index db999028..9cbe18ab 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -53,6 +53,7 @@ DataState &DataState::operator=(DataState &&ss) fix_elem_orient = ss.fix_elem_orient; save_coloring = ss.save_coloring; keep_attr = ss.keep_attr; + cmplx_phase = ss.cmplx_phase; return *this; } @@ -508,11 +509,10 @@ void DataState::SetQuadFunctionSolution(int qf_component) SetQuadSolution(); } -void DataState::SetComplexSolution(ComplexSolution cmplx_type) +void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) { - double (*cmplx_2_scalar)(double, double) = NULL; + double (*cmplx_2_scalar)(double, double); const char *str_cmplx_2_scalar; - GridFunction *gf = NULL; switch (cmplx_type) { @@ -525,11 +525,11 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type) str_cmplx_2_scalar = "phase"; break; case ComplexSolution::Real: - gf = new GridFunction(cgrid_f->FESpace(), cgrid_f->real(), 0); + cmplx_2_scalar = [](double r, double i) { return r; }; str_cmplx_2_scalar = "real part"; break; case ComplexSolution::Imag: - gf = new GridFunction(cgrid_f->FESpace(), cgrid_f->imag(), 0); + cmplx_2_scalar = [](double r, double i) { return i; }; str_cmplx_2_scalar = "imaginary part"; break; default: @@ -537,82 +537,96 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type) return; } - cout << "Representing complex function by: " << str_cmplx_2_scalar << endl; - if (!gf) + if (print) { - gf = new GridFunction(cgrid_f->FESpace()); - const FiniteElementSpace *fes = cgrid_f->FESpace(); - const Mesh *msh = fes->GetMesh(); - const int dim = msh->Dimension(); - if (cgrid_f->FESpace()->FEColl()->GetMapType(dim) == FiniteElement::VALUE) + cout << "Representing complex function by: " << str_cmplx_2_scalar << endl; + } + GridFunction *gf = new GridFunction(cgrid_f->FESpace()); + const FiniteElementSpace *fes = cgrid_f->FESpace(); + const Mesh *msh = fes->GetMesh(); + const int dim = msh->Dimension(); + + 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; + }; + + if (cgrid_f->FESpace()->FEColl()->GetMapType(dim) == FiniteElement::VALUE) + { + for (int i = 0; i < gf->Size(); i++) { - for (int i = 0; i < gf->Size(); i++) - { - (*gf)(i) = cmplx_2_scalar(cgrid_f->real()(i), cgrid_f->imag()(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 + } + else + { + const int vdim = fes->GetVDim(); + const int sdim = fes->GetVectorDim(); + Array vdofs; + ElementTransformation *Tr; + Vector r_data, i_data, z_data; + DenseMatrix r_vals, i_vals, z_vals; + Vector r_vec, i_vec, z_vec; + DenseMatrix vshape; + Vector shape; + for (int z = 0; z < msh->GetNE(); z++) { - const int vdim = fes->GetVDim(); - const int sdim = fes->GetVectorDim(); - Array vdofs; - ElementTransformation *Tr; - Vector r_data, i_data, z_data; - DenseMatrix r_vals, i_vals, z_vals; - Vector r_vec, i_vec, z_vec; - DenseMatrix vshape; - 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); + z_data.SetSize(vdofs.Size()); + + Tr = fes->GetElementTransformation(z); + const FiniteElement *fe = fes->GetFE(z); + const IntegrationRule &ir = fe->GetNodes(); + const int nnp = ir.GetNPoints(); + r_vals.Reset(r_data.GetData(), nnp, vdim); + i_vals.Reset(i_data.GetData(), nnp, vdim); + z_vals.Reset(z_data.GetData(), nnp, vdim); + + if (fe->GetRangeType() == FiniteElement::SCALAR) { - fes->GetElementVDofs(z, vdofs); - cgrid_f->real().GetSubVector(vdofs, r_data); - cgrid_f->imag().GetSubVector(vdofs, i_data); - z_data.SetSize(vdofs.Size()); - - Tr = fes->GetElementTransformation(z); - const FiniteElement *fe = fes->GetFE(z); - const IntegrationRule &ir = fe->GetNodes(); - const int nnp = ir.GetNPoints(); - r_vals.Reset(r_data.GetData(), nnp, vdim); - i_vals.Reset(i_data.GetData(), nnp, vdim); - z_vals.Reset(z_data.GetData(), nnp, vdim); - + shape.SetSize(nnp); + } + else + { + vshape.SetSize(nnp, sdim); + } + for (int n = 0; n < nnp; n++) + { + const IntegrationPoint &ip = ir.IntPoint(n); + Tr->SetIntPoint(&ip); + double w; if (fe->GetRangeType() == FiniteElement::SCALAR) { - shape.SetSize(nnp); + fe->CalcPhysShape(*Tr, shape); + w = shape(n); } else { - vshape.SetSize(nnp, sdim); + fe->CalcPhysVShape(*Tr, vshape); + Vector vec(sdim); + vshape.GetRow(n, vec); + w = vec.Norml2(); } - for (int n = 0; n < nnp; n++) + for (int d = 0; d < vdim; d++) { - const IntegrationPoint &ip = ir.IntPoint(n); - Tr->SetIntPoint(&ip); - double w; - if (fe->GetRangeType() == FiniteElement::SCALAR) - { - fe->CalcPhysShape(*Tr, shape); - w = shape(n); - } - else - { - fe->CalcPhysVShape(*Tr, vshape); - Vector vec(sdim); - vshape.GetRow(n, vec); - w = vec.Norml2(); - } - for (int d = 0; d < vdim; d++) - { - const double rval = r_vals(n,d) * w; - const double ival = i_vals(n,d) * w; - const double zval = cmplx_2_scalar(rval, ival); - z_vals(n,d) = (w != 0.)?(zval / w):(0.); - } + double rval = r_vals(n,d) * w; + double ival = i_vals(n,d) * w; + rot_ph(rval, ival); + const double zval = cmplx_2_scalar(rval, ival); + z_vals(n,d) = (w != 0.)?(zval / w):(0.); } - - gf->SetSubVector(vdofs, z_data); } + + gf->SetSubVector(vdofs, z_data); } } internal.grid_f.reset(gf); @@ -620,9 +634,9 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type) cmplx_sol = cmplx_type; } -void DataState::SwitchComplexSolution(ComplexSolution cmplx_type) +void DataState::SwitchComplexSolution(ComplexSolution cmplx_type, bool print) { - SetComplexSolution(cmplx_type); + SetComplexSolution(cmplx_type, print); ExtrudeMeshAndSolution(); } diff --git a/lib/data_state.hpp b/lib/data_state.hpp index 20ac6efa..7dccd7a2 100644 --- a/lib/data_state.hpp +++ b/lib/data_state.hpp @@ -91,6 +91,7 @@ struct DataState bool fix_elem_orient{false}; bool save_coloring{false}; bool keep_attr{false}; + double cmplx_phase{}; DataState() = default; DataState(DataState &&ss) { *this = std::move(ss); } @@ -176,10 +177,11 @@ struct DataState void SetMeshSolution(); /// Set the complex function representation producing a proxy grid function - void SetComplexSolution(ComplexSolution type = ComplexSolution::Magnitude); + void SetComplexSolution(ComplexSolution type = ComplexSolution::Magnitude, + bool print = true); /// Switch the complex function representation - void SwitchComplexSolution(ComplexSolution type); + void SwitchComplexSolution(ComplexSolution type, bool print = true); /// Get the current representation of complex solution inline ComplexSolution GetComplexSolution() const { return cmplx_sol; } diff --git a/lib/openglvis.hpp b/lib/openglvis.hpp index 5c371063..afa92ec2 100644 --- a/lib/openglvis.hpp +++ b/lib/openglvis.hpp @@ -211,6 +211,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/vsdata.hpp b/lib/vsdata.hpp index 76ab4cd1..b67650fc 100644 --- a/lib/vsdata.hpp +++ b/lib/vsdata.hpp @@ -223,6 +223,8 @@ class VisualizationSceneScalarData : public VisualizationScene virtual gl3::SceneInfo GetSceneObjs(); + virtual void UpdateComplexPhase(double ph) { win.UpdateComplexPhase(ph); } + void ProcessUpdatedBufs(gl3::SceneInfo& scene); void glTF_ExportBox(glTF_Builder &bld, diff --git a/lib/window.cpp b/lib/window.cpp index e7c7ce6e..36422ccf 100644 --- a/lib/window.cpp +++ b/lib/window.cpp @@ -210,6 +210,13 @@ void Window::SwitchQuadSolution(DataState::QuadSolution quad_type) ResetMeshAndSolution(data_state); } +void Window::UpdateComplexPhase(double ph) +{ + data_state.cmplx_phase += ph; + data_state.SwitchComplexSolution(data_state.GetComplexSolution(), false); + ResetMeshAndSolution(data_state); +} + bool Window::SetNewMeshAndSolution(DataState new_state) { if (new_state.mesh->SpaceDimension() == data_state.mesh->SpaceDimension() && diff --git a/lib/window.hpp b/lib/window.hpp index b742e5aa..5cab7228 100644 --- a/lib/window.hpp +++ b/lib/window.hpp @@ -67,6 +67,9 @@ struct Window /// 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. /// From 0dad8ba92716317cb9280e81eaf51a32f40d0c08 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 14 Jan 2026 22:07:04 -0800 Subject: [PATCH 11/44] Added optimization for magnitude. --- lib/window.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/lib/window.cpp b/lib/window.cpp index 36422ccf..09da9d77 100644 --- a/lib/window.cpp +++ b/lib/window.cpp @@ -213,7 +213,10 @@ void Window::SwitchQuadSolution(DataState::QuadSolution quad_type) void Window::UpdateComplexPhase(double ph) { data_state.cmplx_phase += ph; - data_state.SwitchComplexSolution(data_state.GetComplexSolution(), false); + DataState::ComplexSolution cs = data_state.GetComplexSolution(); + // check if magnitude is viewed, which remains the same + if (cs == DataState::ComplexSolution::Magnitude) { return; } + data_state.SwitchComplexSolution(cs, false); ResetMeshAndSolution(data_state); } From cde5753aeacabf480d5e415a02e91f4c01396090 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 09:36:30 -0800 Subject: [PATCH 12/44] Minor clarification of phase rate. --- lib/aux_vis.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/aux_vis.cpp b/lib/aux_vis.cpp index 3380441a..8e6dc51e 100644 --- a/lib/aux_vis.cpp +++ b/lib/aux_vis.cpp @@ -1464,7 +1464,7 @@ void CheckPhaseAnim() phase_anim = 0; RemoveIdleFunc(MainLoop); } - cout << "Phase rate: " << phase_rate << " radians / frame" << endl; + cout << "Phase rate: " << phase_rate << " period / frame" << endl; } void KeySemicolonPressed() From 67ad9ebf0c2cff97a3c2483e6ee884ceb75077bb Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 09:40:33 -0800 Subject: [PATCH 13/44] Prevented reduction of precision in phase. --- lib/window.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/window.cpp b/lib/window.cpp index 09da9d77..54d6b900 100644 --- a/lib/window.cpp +++ b/lib/window.cpp @@ -213,6 +213,7 @@ void Window::SwitchQuadSolution(DataState::QuadSolution quad_type) 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; } From c709ba92ce762d0711d39f84d3440564d25c8692 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 09:58:15 -0800 Subject: [PATCH 14/44] Decoupled phase animation from spinning. --- lib/aux_vis.cpp | 42 +++++++++++++++++++----------------------- lib/aux_vis.hpp | 1 + 2 files changed, 20 insertions(+), 23 deletions(-) diff --git a/lib/aux_vis.cpp b/lib/aux_vis.cpp index 8e6dc51e..11351225 100644 --- a/lib/aux_vis.cpp +++ b/lib/aux_vis.cpp @@ -619,6 +619,18 @@ thread_local double phase_rate = 0.; thread_local int phase_anim = 0; const double phase_step = 0.01; +void CheckMainIdleFunc() +{ + if (locscene->spinning || phase_anim) + { + AddIdleFunc(MainLoop); + } + if (!locscene->spinning && !phase_anim) + { + RemoveIdleFunc(MainLoop); + } +} + void MainLoop() { static int p = 1; @@ -683,7 +695,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; @@ -1395,16 +1407,8 @@ 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; } @@ -1426,14 +1430,14 @@ void KeyDeletePressed() { xang = yang = 0.; locscene -> spinning = 0; - RemoveIdleFunc(MainLoop); + CheckMainIdleFunc(); constrained_spinning = 1; } else { xang = xang_step; locscene -> spinning = 1; - AddIdleFunc(MainLoop); + CheckMainIdleFunc(); constrained_spinning = 1; } } @@ -1454,16 +1458,8 @@ void CheckPhaseAnim() { phase_rate = 0.; } - if (phase_rate != 0.) - { - phase_anim = 1; - AddIdleFunc(MainLoop); - } - else - { - phase_anim = 0; - RemoveIdleFunc(MainLoop); - } + phase_anim = (phase_rate != 0.); + CheckMainIdleFunc(); cout << "Phase rate: " << phase_rate << " period / frame" << endl; } diff --git a/lib/aux_vis.hpp b/lib/aux_vis.hpp index c77fc2b7..c9d40caa 100644 --- a/lib/aux_vis.hpp +++ b/lib/aux_vis.hpp @@ -53,6 +53,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); From 638921e4ad681c05e4f9d0ab319abf52347e2e8d Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 10:49:28 -0800 Subject: [PATCH 15/44] Replaced the stream command for complex solutions by automatic detection. --- lib/file_reader.cpp | 2 +- lib/file_reader.hpp | 2 +- lib/stream_reader.cpp | 53 +++++++++++++++++++++++++++++-------------- lib/stream_reader.hpp | 2 ++ 4 files changed, 40 insertions(+), 19 deletions(-) diff --git a/lib/file_reader.cpp b/lib/file_reader.cpp index 27d5a23c..a2662b84 100644 --- a/lib/file_reader.cpp +++ b/lib/file_reader.cpp @@ -128,7 +128,7 @@ int FileReader::ReadParallel(int np, FileType ft, const char *mesh_file, return read_err; } -bool FileReader::CheckStreamIsComplex(std::istream &solin, bool parallel) const +bool FileReader::CheckStreamIsComplex(std::istream &solin, bool parallel) { string buff; auto pos = solin.tellg(); diff --git a/lib/file_reader.hpp b/lib/file_reader.hpp index 337704cd..5a58b484 100644 --- a/lib/file_reader.hpp +++ b/lib/file_reader.hpp @@ -21,7 +21,7 @@ class FileReader DataState &data; int pad_digits; - bool CheckStreamIsComplex(std::istream &sol, bool parallel = false) const; + static bool CheckStreamIsComplex(std::istream &sol, bool parallel = false); int ReadParMeshAndGridFunction(int np, const char *mesh_prefix, const char *sol_prefix, int component = -1); diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 5d5196bb..3564f1d8 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -13,6 +13,7 @@ #include #include +#include using namespace std; using namespace mfem; @@ -21,7 +22,6 @@ enum class Command { Mesh, Solution, - ComplexSolution, Quadrature, Fem2D, VFem2D, @@ -69,7 +69,6 @@ StreamCommands::StreamCommands() { (*this)[Command::Mesh] = {"mesh", false, "", "Visualize the mesh."}; (*this)[Command::Solution] = {"solution", false, " ", "Visualize the solution."}; - (*this)[Command::ComplexSolution] = {"csolution", false, " ", "Visualize the complex solution."}; (*this)[Command::Quadrature] = {"quadrature", false, " ", "Visualize the quadrature."}; (*this)[Command::Fem2D] = {"fem2d_data", false, " ", "Visualize the 2D scalar data."}; (*this)[Command::VFem2D] = {"vfem2d_data", false, " ", "Visualize the 2D vector data."}; @@ -104,6 +103,21 @@ bool StreamReader::SupportsDataType(const string &data_type) return it != commands.end(); } +bool StreamReader::CheckStreamIsComplex(std::istream &solin, + bool parallel) +{ + string buff; + solin >> ws; + getline(solin, buff); + solin.unget(); + const size_t len = buff.length(); + for (size_t i = 0; i < len; i++) { solin.putback(buff[len-1-i]); } + mfem::filter_dos(buff); + const char *header = (parallel)?("ParComplexGridFunction"): + ("ComplexGridFunction"); + return (buff == header); +} + int StreamReader::ReadStream( istream &is, const string &data_type) { @@ -156,11 +170,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)); - break; - case Command::ComplexSolution: - data.SetMesh(new Mesh(is, 1, 0, data.fix_elem_orient)); - data.SetCmplxGridFunction(new ComplexGridFunction(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)); @@ -338,22 +355,24 @@ int StreamReader::ReadStreams(const StreamCollection &input_streams) } } gf_array[p] = NULL; - if (data_type == "csolution") - { - isock >> ws; - isock.ignore(3);// ignore 'Par' prefix to load as serial - cgf_array[p] = new ComplexGridFunction(mesh_array[p], isock); - cgf_count++; - } - else if (data_type == "quadrature") + if (data_type == "quadrature") { qf_array[p] = new QuadratureFunction(mesh_array[p], isock); qf_count++; } 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; diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index 5801d362..f9eb2d67 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -25,6 +25,8 @@ class StreamReader { DataState &data; + static bool CheckStreamIsComplex(std::istream &sol, bool parallel = false); + public: StreamReader(DataState &data_): data(data_) { } From 8cc445e28fa8b3e9708a4b46f28649aa856ff91a Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 10:56:12 -0800 Subject: [PATCH 16/44] Reduced the default phase step. --- lib/aux_vis.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/aux_vis.cpp b/lib/aux_vis.cpp index 11351225..f1028c06 100644 --- a/lib/aux_vis.cpp +++ b/lib/aux_vis.cpp @@ -617,7 +617,7 @@ thread_local int constrained_spinning = 0; thread_local double phase_rate = 0.; thread_local int phase_anim = 0; -const double phase_step = 0.01; +const double phase_step = 0.001; void CheckMainIdleFunc() { @@ -1454,7 +1454,7 @@ void KeyEnterPressed() void CheckPhaseAnim() { - if (fabs(phase_rate) < 1.e-2) + if (fabs(phase_rate) < phase_step) { phase_rate = 0.; } From a532240f54362c205ace35ec35e18ec2892baea4 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 11:02:32 -0800 Subject: [PATCH 17/44] Fixed the phas rate check. --- lib/aux_vis.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/aux_vis.cpp b/lib/aux_vis.cpp index f1028c06..19810ea9 100644 --- a/lib/aux_vis.cpp +++ b/lib/aux_vis.cpp @@ -1454,7 +1454,7 @@ void KeyEnterPressed() void CheckPhaseAnim() { - if (fabs(phase_rate) < phase_step) + if (fabs(phase_rate) < phase_step / 2.) { phase_rate = 0.; } From 5cd82611f7d2bfb648754f4773039b1042a9422f Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 12:53:19 -0800 Subject: [PATCH 18/44] Updated Changelog. --- CHANGELOG | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/CHANGELOG b/CHANGELOG index fd116769..397547f9 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. + Version 4.4 released on May 1, 2025 =================================== From 53ca05dc686e3672936c1c63e346ad6c7f40891a Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 14:09:37 -0800 Subject: [PATCH 19/44] Added documentation of the keys with Q for cycling for now. --- README.md | 6 ++++++ lib/vssolution.cpp | 2 ++ lib/vssolution3d.cpp | 2 ++ lib/vsvector.cpp | 2 ++ lib/vsvector3d.cpp | 2 ++ 5 files changed, 14 insertions(+) diff --git a/README.md b/README.md index ee4c3dee..18c5a572 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) +- Q – Cycle between representations of the visualized *complex data*. The options are: + - magnitude + - phase + - real part + - imaginary part - q – Exit ## Advanced @@ -199,6 +204,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 +- ; / : - Increase/decrease the harmonic complex phase rate of animation ## 2D scalar data diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 44638f00..604b747b 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -76,6 +76,8 @@ std::string VisualizationSceneSolution::GetHelpString() const << "| p/P Cycle through color palettes |" << endl << "| q - Quits |" << endl << "| Q - Cycle quadrature data mode |" << endl + << "| Q - Cycle complex data mode |" << endl + << "| ;/: - Inc/decrease phase freq. |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index ada41aec..3337e4cc 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -65,6 +65,8 @@ std::string VisualizationSceneSolution3d::GetHelpString() const << "| p/P Cycle through color palettes |" << endl << "| q - Quits |" << endl << "| Q - Cycle quadrature data mode |" << endl + << "| Q - Cycle complex data mode |" << endl + << "| ;/: - Inc/decrease phase freq. |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index 206f47fd..84139c47 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -52,6 +52,8 @@ std::string VisualizationSceneVector::GetHelpString() const << "| p/P Cycle through color palettes |" << endl << "| q - Quits |" << endl << "| Q - Cycle quadrature data mode |" << endl + << "| Q - Cycle complex data mode |" << endl + << "| ;/: - Inc/decrease phase freq. |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index c9f99f90..03e25af9 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -62,6 +62,8 @@ std::string VisualizationSceneVector3d::GetHelpString() const << "| p/P Cycle through color palettes |" << endl << "| q - Quits |" << endl << "| Q - Cycle quadrature data mode |" << endl + << "| Q - Cycle complex data mode |" << endl + << "| ;/: - Inc/decrease phase freq. |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl From c9612739abe427a173de55337338e4ea60d1d671 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 16:34:26 -0800 Subject: [PATCH 20/44] Fixed the value ranges for phase animation. --- lib/data_state.cpp | 6 ++++++ lib/data_state.hpp | 1 + lib/vsdata.cpp | 7 +++++-- lib/vsdata.hpp | 2 +- lib/window.cpp | 20 ++++++++++++++++++++ 5 files changed, 33 insertions(+), 3 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index 9cbe18ab..0f5709d0 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -555,12 +555,16 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) r = rph, i = iph; }; + cmplx_mag_max = 0.; + if (cgrid_f->FESpace()->FEColl()->GetMapType(dim) == FiniteElement::VALUE) { for (int i = 0; i < gf->Size(); i++) { double rval = cgrid_f->real()(i); double ival = cgrid_f->imag()(i); + double mag = hypot(rval, ival); + cmplx_mag_max = max(cmplx_mag_max, mag); rot_ph(rval, ival); (*gf)(i) = cmplx_2_scalar(rval, ival); } @@ -620,6 +624,8 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) { double rval = r_vals(n,d) * w; double ival = i_vals(n,d) * w; + double mag = hypot(rval, ival); + cmplx_mag_max = min(cmplx_mag_max, mag); rot_ph(rval, ival); const double zval = cmplx_2_scalar(rval, ival); z_vals(n,d) = (w != 0.)?(zval / w):(0.); diff --git a/lib/data_state.hpp b/lib/data_state.hpp index 7dccd7a2..a633f452 100644 --- a/lib/data_state.hpp +++ b/lib/data_state.hpp @@ -92,6 +92,7 @@ struct DataState bool save_coloring{false}; bool keep_attr{false}; double cmplx_phase{}; + double cmplx_mag_max{0.}; DataState() = default; DataState(DataState &&ss) { *this = std::move(ss); } diff --git a/lib/vsdata.cpp b/lib/vsdata.cpp index f662b14a..24adcc2b 100644 --- a/lib/vsdata.cpp +++ b/lib/vsdata.cpp @@ -1346,12 +1346,15 @@ void VisualizationSceneScalarData::ToggleTexture() } } -void VisualizationSceneScalarData::SetAutoscale(int _autoscale) +void VisualizationSceneScalarData::SetAutoscale(int _autoscale, bool update) { if (autoscale != _autoscale) { autoscale = _autoscale; - DoAutoscale(true); + if (update) + { + DoAutoscale(true); + } } } diff --git a/lib/vsdata.hpp b/lib/vsdata.hpp index b67650fc..79e6a0bc 100644 --- a/lib/vsdata.hpp +++ b/lib/vsdata.hpp @@ -316,7 +316,7 @@ class VisualizationSceneScalarData : public VisualizationScene void Toggle2DView(); - void SetAutoscale(int _autoscale); + void SetAutoscale(int _autoscale, bool update = true); int GetAutoscale() const { return autoscale; } /// Shrink the set of points towards attributes centers of gravity diff --git a/lib/window.cpp b/lib/window.cpp index 54d6b900..998720fb 100644 --- a/lib/window.cpp +++ b/lib/window.cpp @@ -201,7 +201,23 @@ void Window::GLVisStartVis() void Window::SwitchComplexSolution(DataState::ComplexSolution cmplx_type) { data_state.SwitchComplexSolution(cmplx_type); + int as = (cmplx_type == DataState::ComplexSolution::Magnitude ? 1 : 3); + vs->SetAutoscale(as, false); ResetMeshAndSolution(data_state); + switch (cmplx_type) + { + case DataState::ComplexSolution::Magnitude: + break; + case DataState::ComplexSolution::Phase: + vs->SetValueRange(-M_PI, +M_PI); + break; + case DataState::ComplexSolution::Real: + case DataState::ComplexSolution::Imag: + vs->SetValueRange(-data_state.cmplx_mag_max, +data_state.cmplx_mag_max); + break; + default: + mfem_error("Unknown representation"); + } } void Window::SwitchQuadSolution(DataState::QuadSolution quad_type) @@ -218,7 +234,11 @@ void Window::UpdateComplexPhase(double ph) // check if magnitude is viewed, which remains the same if (cs == DataState::ComplexSolution::Magnitude) { return; } data_state.SwitchComplexSolution(cs, false); + // do not autoscale for animation + int as = vs->GetAutoscale(); + vs->SetAutoscale(0); ResetMeshAndSolution(data_state); + vs->SetAutoscale(as, false); } bool Window::SetNewMeshAndSolution(DataState new_state) From cd2b2daca71702760d5076902464355f3811792d Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 16:38:31 -0800 Subject: [PATCH 21/44] Added a note about animated phase. --- lib/data_state.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index 0f5709d0..0f50a6b4 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -539,7 +539,9 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) if (print) { - cout << "Representing complex function by: " << str_cmplx_2_scalar << endl; + cout << "Representing complex function by: " << str_cmplx_2_scalar; + if (cmplx_type != ComplexSolution::Magnitude) { cout << " (+ animated harmonic phase)"; } + cout << endl; } GridFunction *gf = new GridFunction(cgrid_f->FESpace()); const FiniteElementSpace *fes = cgrid_f->FESpace(); From d7d4edd93b5b1de1943d1f807377faa03ac653b9 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 15 Jan 2026 18:26:16 -0800 Subject: [PATCH 22/44] Fixed a typo in max of complex magnitude. --- lib/data_state.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index 0f50a6b4..8a6c2c70 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -627,7 +627,7 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) double rval = r_vals(n,d) * w; double ival = i_vals(n,d) * w; double mag = hypot(rval, ival); - cmplx_mag_max = min(cmplx_mag_max, mag); + cmplx_mag_max = max(cmplx_mag_max, mag); rot_ph(rval, ival); const double zval = cmplx_2_scalar(rval, ival); z_vals(n,d) = (w != 0.)?(zval / w):(0.); From 4dad4dd145edad82976fe530c047710d3d5e6196 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 16 Jan 2026 11:02:04 -0800 Subject: [PATCH 23/44] Refactored autoscale modes. --- lib/script_controller.cpp | 8 ++++---- lib/threads.cpp | 8 ++++---- lib/vsdata.cpp | 22 ++++++++++++---------- lib/vsdata.hpp | 33 +++++++++++++++++++++------------ lib/vsvector.cpp | 2 +- lib/window.cpp | 10 ++++++---- 6 files changed, 48 insertions(+), 35 deletions(-) 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/threads.cpp b/lib/threads.cpp index a0aa6509..d5cfe56d 100644 --- a/lib/threads.cpp +++ b/lib/threads.cpp @@ -694,19 +694,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 24adcc2b..8eece3ff 100644 --- a/lib/vsdata.cpp +++ b/lib/vsdata.cpp @@ -105,15 +105,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); } @@ -121,7 +121,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); } @@ -666,11 +666,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(); } @@ -1346,7 +1347,8 @@ void VisualizationSceneScalarData::ToggleTexture() } } -void VisualizationSceneScalarData::SetAutoscale(int _autoscale, bool update) +void VisualizationSceneScalarData::SetAutoscale(Autoscale _autoscale, + bool update) { if (autoscale != _autoscale) { @@ -1473,7 +1475,7 @@ void VisualizationSceneScalarData::Init() PrepareRuler(); - autoscale = 1; + autoscale = VisualizationSceneScalarData::Autoscale::MeshAndValue; } VisualizationSceneScalarData::~VisualizationSceneScalarData() diff --git a/lib/vsdata.hpp b/lib/vsdata.hpp index 79e6a0bc..729a2336 100644 --- a/lib/vsdata.hpp +++ b/lib/vsdata.hpp @@ -66,6 +66,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: Mesh *mesh{}, *mesh_coarse{}; Vector *sol{}; @@ -104,12 +118,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; @@ -187,10 +196,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); @@ -316,8 +325,8 @@ class VisualizationSceneScalarData : public VisualizationScene void Toggle2DView(); - void SetAutoscale(int _autoscale, bool update = true); - 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(DenseMatrix &pointmat, int i, int fn, int di); diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index 84139c47..c2a957a2 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -459,7 +459,7 @@ void VisualizationSceneVector::NewMeshAndSolution(GridFunction &vgf, Mesh *mc) VisualizationSceneSolution::NewMeshAndSolution(mesh, mesh_coarse, sol, &vgf); - if (autoscale) + if (autoscale != VisualizationSceneScalarData::Autoscale::None) { if (Vec2Scalar == VecLength) { diff --git a/lib/window.cpp b/lib/window.cpp index 998720fb..893aebfa 100644 --- a/lib/window.cpp +++ b/lib/window.cpp @@ -167,7 +167,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) @@ -201,7 +201,9 @@ void Window::GLVisStartVis() void Window::SwitchComplexSolution(DataState::ComplexSolution cmplx_type) { data_state.SwitchComplexSolution(cmplx_type); - int as = (cmplx_type == DataState::ComplexSolution::Magnitude ? 1 : 3); + auto as = (cmplx_type == DataState::ComplexSolution::Magnitude ? + VisualizationSceneScalarData::Autoscale::MeshAndValue : + VisualizationSceneScalarData::Autoscale::Mesh); vs->SetAutoscale(as, false); ResetMeshAndSolution(data_state); switch (cmplx_type) @@ -235,8 +237,8 @@ void Window::UpdateComplexPhase(double ph) if (cs == DataState::ComplexSolution::Magnitude) { return; } data_state.SwitchComplexSolution(cs, false); // do not autoscale for animation - int as = vs->GetAutoscale(); - vs->SetAutoscale(0); + auto as = vs->GetAutoscale(); + vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::None); ResetMeshAndSolution(data_state); vs->SetAutoscale(as, false); } From 60d571191b165324efa94bc2d9d400c031894244 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 23 Jan 2026 13:19:51 -0800 Subject: [PATCH 24/44] Changed phase animation keys to Alt+0/./Enter. --- README.md | 2 +- lib/aux_vis.cpp | 106 ++++++++++++++++++++++++------------------- lib/aux_vis.hpp | 9 ++-- lib/vssolution.cpp | 3 +- lib/vssolution3d.cpp | 3 +- lib/vsvector.cpp | 3 +- lib/vsvector3d.cpp | 3 +- 7 files changed, 72 insertions(+), 57 deletions(-) diff --git a/README.md b/README.md index 18c5a572..9e417378 100644 --- a/README.md +++ b/README.md @@ -204,7 +204,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 -- ; / : - Increase/decrease the harmonic complex phase rate of animation +- Alt+. – Start/stop harmonic complex phase animation (speed can be controlled with Alt+0 / Alt+Enter) ## 2D scalar data diff --git a/lib/aux_vis.cpp b/lib/aux_vis.cpp index 19810ea9..c2c2403d 100644 --- a/lib/aux_vis.cpp +++ b/lib/aux_vis.cpp @@ -229,9 +229,6 @@ GLWindow* InitVisualization(const char name[], int x, int y, int w, int h, wnd->setOnKeyDown (SDLK_PERIOD, KeyDeletePressed); wnd->setOnKeyDown (SDLK_RETURN, KeyEnterPressed); - wnd->setOnKeyDown (SDLK_SEMICOLON, KeySemicolonPressed); - wnd->setOnKeyDown (SDLK_COLON, KeyColonPressed); - wnd->setOnKeyDown (SDLK_0, Key0Pressed); wnd->setOnKeyDown (SDLK_1, Key1Pressed); wnd->setOnKeyDown (SDLK_2, Key2Pressed); @@ -609,6 +606,7 @@ 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; @@ -1412,75 +1410,91 @@ void CheckSpin() cout << "Spin angle: " << xang << " degrees / frame" << endl; } -const double xang_step = 0.2; // angle in degrees - -void Key0Pressed() +void CheckPhaseAnim() { - if (!locscene -> spinning) + if (fabs(phase_rate) < phase_step / 2.) { - xang = 0; + phase_rate = 0.; } - xang -= xang_step; - CheckSpin(); + phase_anim = (phase_rate != 0.); + CheckMainIdleFunc(); + cout << "Phase rate: " << phase_rate << " period / frame" << endl; } -void KeyDeletePressed() +void Key0Pressed(SDL_Keymod mod) { - if (locscene -> spinning) + if (mod & KMOD_ALT) { - xang = yang = 0.; - locscene -> spinning = 0; - CheckMainIdleFunc(); - constrained_spinning = 1; + if (!phase_anim) + { + phase_rate = 0.; + } + phase_rate -= phase_step; + CheckPhaseAnim(); } else { - xang = xang_step; - locscene -> spinning = 1; - CheckMainIdleFunc(); - constrained_spinning = 1; + if (!locscene -> spinning) + { + xang = 0; + } + xang -= xang_step; + CheckSpin(); } } -void KeyEnterPressed() +void KeyDeletePressed(SDL_Keymod mod) { - if (!locscene -> spinning) + if (mod & KMOD_ALT) { - xang = 0; + if (phase_anim) + { + phase_rate = 0.; + phase_anim = 0; + } + else + { + phase_rate = phase_step; + phase_anim = 1; + } } - xang += xang_step; - CheckSpin(); -} - -void CheckPhaseAnim() -{ - if (fabs(phase_rate) < phase_step / 2.) + else { - phase_rate = 0.; + if (locscene -> spinning) + { + xang = yang = 0.; + locscene -> spinning = 0; + } + else + { + xang = xang_step; + locscene -> spinning = 1; + } + constrained_spinning = 1; } - phase_anim = (phase_rate != 0.); CheckMainIdleFunc(); - cout << "Phase rate: " << phase_rate << " period / frame" << endl; } -void KeySemicolonPressed() +void KeyEnterPressed(SDL_Keymod mod) { - if (!phase_anim) + if (mod & KMOD_ALT) { - phase_rate = 0.; + if (!phase_anim) + { + phase_rate = 0.; + } + phase_rate += phase_step; + CheckPhaseAnim(); } - phase_rate += phase_step; - CheckPhaseAnim(); -} - -void KeyColonPressed() -{ - if (!phase_anim) + else { - phase_rate = 0.; + if (!locscene -> spinning) + { + xang = 0; + } + xang += xang_step; + CheckSpin(); } - phase_rate -= phase_step; - CheckPhaseAnim(); } void Key7Pressed() diff --git a/lib/aux_vis.hpp b/lib/aux_vis.hpp index c9d40caa..1b96ac1a 100644 --- a/lib/aux_vis.hpp +++ b/lib/aux_vis.hpp @@ -85,12 +85,9 @@ void Key7Pressed(); void Key8Pressed(); void Key9Pressed(); -void Key0Pressed(); -void KeyDeletePressed(); -void KeyEnterPressed(); - -void KeySemicolonPressed(); -void KeyColonPressed(); +void Key0Pressed(SDL_Keymod); +void KeyDeletePressed(SDL_Keymod); +void KeyEnterPressed(SDL_Keymod); void KeyLeftPressed(GLenum); void KeyRightPressed(GLenum); diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 604b747b..247c7679 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -77,7 +77,6 @@ std::string VisualizationSceneSolution::GetHelpString() const << "| q - Quits |" << endl << "| Q - Cycle quadrature data mode |" << endl << "| Q - Cycle complex data mode |" << endl - << "| ;/: - Inc/decrease phase freq. |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl @@ -110,6 +109,8 @@ std::string VisualizationSceneSolution::GetHelpString() const << "| +/- Change z-scaling |" << endl << "| . - Start/stop spinning |" << endl << "| 0/Enter - Spinning speed and dir. |" << endl + << "| Alt+. - Start/stop phase anim. |" << endl + << "| Alt+0/Enter - Phase anim. speed |" << endl << "+------------------------------------+" << endl << "| Mouse |" << endl << "+------------------------------------+" << endl diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 3337e4cc..0293774f 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -66,7 +66,6 @@ std::string VisualizationSceneSolution3d::GetHelpString() const << "| q - Quits |" << endl << "| Q - Cycle quadrature data mode |" << endl << "| Q - Cycle complex data mode |" << endl - << "| ;/: - Inc/decrease phase freq. |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl @@ -104,6 +103,8 @@ std::string VisualizationSceneSolution3d::GetHelpString() const << "| +/- Change z-scaling |" << endl << "| . - Start/stop spinning |" << endl << "| 0/Enter - Spinning speed and dir. |" << endl + << "| Alt+. - Start/stop phase anim. |" << endl + << "| Alt+0/Enter - Phase anim. speed |" << endl << "+------------------------------------+" << endl << "| Mouse |" << endl << "+------------------------------------+" << endl diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index c2a957a2..e78364a6 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -53,7 +53,6 @@ std::string VisualizationSceneVector::GetHelpString() const << "| q - Quits |" << endl << "| Q - Cycle quadrature data mode |" << endl << "| Q - Cycle complex data mode |" << endl - << "| ;/: - Inc/decrease phase freq. |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl @@ -87,6 +86,8 @@ std::string VisualizationSceneVector::GetHelpString() const << "| +/- Change z-scaling |" << endl << "| . - Start/stop spinning |" << endl << "| 0/Enter - Spinning speed and dir. |" << endl + << "| Alt+. - Start/stop phase anim. |" << endl + << "| Alt+0/Enter - Phase anim. speed |" << endl << "+------------------------------------+" << endl << "| Mouse |" << endl << "+------------------------------------+" << endl diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index 03e25af9..fbbdb36d 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -63,7 +63,6 @@ std::string VisualizationSceneVector3d::GetHelpString() const << "| q - Quits |" << endl << "| Q - Cycle quadrature data mode |" << endl << "| Q - Cycle complex data mode |" << endl - << "| ;/: - Inc/decrease phase freq. |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl @@ -99,6 +98,8 @@ std::string VisualizationSceneVector3d::GetHelpString() const << "| +/- Change z-scaling |" << endl << "| . - Start/stop spinning |" << endl << "| 0/Enter - Spinning speed and dir. |" << endl + << "| Alt+. - Start/stop phase anim. |" << endl + << "| Alt+0/Enter - Phase anim. speed |" << endl << "+------------------------------------+" << endl << "| Mouse |" << endl << "+------------------------------------+" << endl From 745f96eb6b549a95a58b12e8bb17c51e963bca88 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 23 Jan 2026 13:22:41 -0800 Subject: [PATCH 25/44] Minor cleanup --- lib/aux_vis.cpp | 10 +++++----- lib/aux_vis.hpp | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/lib/aux_vis.cpp b/lib/aux_vis.cpp index c2c2403d..7204ff70 100644 --- a/lib/aux_vis.cpp +++ b/lib/aux_vis.cpp @@ -1369,7 +1369,7 @@ void ToggleThreads() } } -void ThreadsPauseFunc(GLenum state) +void ThreadsPauseFunc(SDL_Keymod state) { if (state & KMOD_CTRL) { @@ -1573,7 +1573,7 @@ void ShiftView(double dx, double dy) locscene->ViewCenterY += dy/scale; } -void KeyLeftPressed(GLenum state) +void KeyLeftPressed(SDL_Keymod state) { if (state & KMOD_CTRL) { @@ -1586,7 +1586,7 @@ void KeyLeftPressed(GLenum state) SendExposeEvent(); } -void KeyRightPressed(GLenum state) +void KeyRightPressed(SDL_Keymod state) { if (state & KMOD_CTRL) { @@ -1599,7 +1599,7 @@ void KeyRightPressed(GLenum state) SendExposeEvent(); } -void KeyUpPressed(GLenum state) +void KeyUpPressed(SDL_Keymod state) { if (state & KMOD_CTRL) { @@ -1612,7 +1612,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 1b96ac1a..8446f85e 100644 --- a/lib/aux_vis.hpp +++ b/lib/aux_vis.hpp @@ -71,7 +71,7 @@ void KeyCtrlP(); void KeyS(); void KeyQPressed(); void ToggleThreads(); -void ThreadsPauseFunc(GLenum); +void ThreadsPauseFunc(SDL_Keymod); void ThreadsStop(); void ThreadsRun(); @@ -89,10 +89,10 @@ 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(); From 59ba959861f4aa3b06d75f46c22f90070495bcae Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 23 Jan 2026 17:19:32 -0800 Subject: [PATCH 26/44] Fixed README for complex phase anim. --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 9e417378..f6f8d197 100644 --- a/README.md +++ b/README.md @@ -204,7 +204,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 -- Alt+. – Start/stop harmonic complex phase animation (speed can be controlled with Alt+0 / Alt+Enter) +- Alt + . – Start/stop harmonic complex phase animation (speed can be controlled with Alt + 0 / Alt + Enter) ## 2D scalar data From 60c38025551cb91174939909d026fb7cd6c484cc Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 23 Jan 2026 19:22:54 -0800 Subject: [PATCH 27/44] Added switching of 3D complex vector grid functions to DGs. --- lib/data_state.cpp | 37 +++++++++++++++++++++++++++++++++++++ lib/data_state.hpp | 6 ++++-- 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index 795952ee..59dc0da4 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -856,6 +856,43 @@ DataState::ProjectVectorFEGridFunction(std::unique_ptr gf) return gf; } +std::unique_ptr +DataState::ProjectVectorFEGridFunction(std::unique_ptr 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)); + SwitchComplexSolution(GetComplexSolution(), false); + } + else + { + internal.grid_f = ProjectVectorFEGridFunction(std::move(internal.grid_f)); + } +} + void ::VectorExtrudeCoefficient::Eval(Vector &v, ElementTransformation &T, const IntegrationPoint &ip) { diff --git a/lib/data_state.hpp b/lib/data_state.hpp index e9b1e933..ee9c257b 100644 --- a/lib/data_state.hpp +++ b/lib/data_state.hpp @@ -218,8 +218,10 @@ struct DataState static std::unique_ptr ProjectVectorFEGridFunction(std::unique_ptr gf); - void ProjectVectorFEGridFunction() - { internal.grid_f = ProjectVectorFEGridFunction(std::move(internal.grid_f)); } + static std::unique_ptr + ProjectVectorFEGridFunction(std::unique_ptr gf); + + void ProjectVectorFEGridFunction(); }; #endif // GLVIS_DATA_STATE_HPP From 4bc003589ca7cb0b702352a06e90c39c6859ee91 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 26 Jan 2026 12:57:24 -0800 Subject: [PATCH 28/44] Generalized offsets for numberings to complex grid functions. --- lib/data_state.cpp | 31 ++++++++++++++++++++++++------- lib/data_state.hpp | 2 +- 2 files changed, 25 insertions(+), 8 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index 353bbde1..d30c1ea9 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -161,7 +161,15 @@ void DataState::SetGridFunction(std::vector &gf_array, { SetGridFunction(new GridFunction(mesh.get(), gf_array.data(), num_pieces), component); - if (!keep_attr) { ComputeDofsOffsets(gf_array); } + if (!keep_attr) + { + std::vector 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, @@ -213,6 +221,15 @@ void DataState::SetCmplxGridFunction( delete rgf; delete igf; SetCmplxGridFunction(cgf, component); + if (!keep_attr) + { + std::vector 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) @@ -902,13 +919,14 @@ void DataState::ProjectVectorFEGridFunction() } } -void DataState::ComputeDofsOffsets(std::vector &gf_array) +void DataState::ComputeDofsOffsets(std::vector + &fespaces) { - const int nprocs = static_cast(gf_array.size()); - MFEM_VERIFY(!gf_array.empty(), "No grid functions provided for offsets"); + const int nprocs = static_cast(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(nprocs); @@ -916,8 +934,7 @@ void DataState::ComputeDofsOffsets(std::vector &gf_array) Array 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 986f6117..4ed8f0d2 100644 --- a/lib/data_state.hpp +++ b/lib/data_state.hpp @@ -104,7 +104,7 @@ struct DataState void SetQuadFunctionSolution(int component = -1); /// Compute the dofs offsets from the grid function vector - void ComputeDofsOffsets(std::vector &gf_array); + void ComputeDofsOffsets(std::vector &fespaces); public: const std::unique_ptr &sol{internal.sol}; From c6d38c693b1cc85f6e001b77697eebcb9073007d Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 26 Jan 2026 13:40:18 -0800 Subject: [PATCH 29/44] Fixed override of UpdateComplexPhase(). --- lib/vsdata.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/vsdata.hpp b/lib/vsdata.hpp index 3bc4bac3..3c55d3e7 100644 --- a/lib/vsdata.hpp +++ b/lib/vsdata.hpp @@ -233,7 +233,7 @@ class VisualizationSceneScalarData : public VisualizationScene gl3::SceneInfo GetSceneObjs() override; - virtual void UpdateComplexPhase(double ph) { win.UpdateComplexPhase(ph); } + void UpdateComplexPhase(double ph) override { win.UpdateComplexPhase(ph); } void ProcessUpdatedBufs(gl3::SceneInfo& scene); From 56eb4cd36aac93a466d5be1223f69dece9f76c8c Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 26 Jan 2026 13:47:58 -0800 Subject: [PATCH 30/44] Fixed CMake compilation. --- lib/file_reader.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/file_reader.cpp b/lib/file_reader.cpp index f52156c2..b4844a01 100644 --- a/lib/file_reader.cpp +++ b/lib/file_reader.cpp @@ -9,11 +9,11 @@ // terms of the BSD-3 license. We welcome feedback and contributions, see file // CONTRIBUTING.md for details. +#include "file_reader.hpp" + #include #include -#include "file_reader.hpp" - using namespace std; using namespace mfem; From 40a72959f7b425bec999ef8e2c20ef5e9f6f96d8 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 26 Jan 2026 17:58:20 -0800 Subject: [PATCH 31/44] Minor DataState cleanup. --- lib/data_state.hpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/lib/data_state.hpp b/lib/data_state.hpp index 4ed8f0d2..a93ab637 100644 --- a/lib/data_state.hpp +++ b/lib/data_state.hpp @@ -103,6 +103,12 @@ struct DataState void SetComplexFunctionSolution(int component = -1); void SetQuadFunctionSolution(int component = -1); + static std::unique_ptr + ProjectVectorFEGridFunction(std::unique_ptr gf); + + static std::unique_ptr + ProjectVectorFEGridFunction(std::unique_ptr gf); + /// Compute the dofs offsets from the grid function vector void ComputeDofsOffsets(std::vector &fespaces); @@ -250,12 +256,6 @@ 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 - ProjectVectorFEGridFunction(std::unique_ptr gf); - - static std::unique_ptr - ProjectVectorFEGridFunction(std::unique_ptr gf); - void ProjectVectorFEGridFunction(); }; From e072e4564068e423f53434807e7634ad0e2ea6fd Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 27 Jan 2026 13:01:35 -0800 Subject: [PATCH 32/44] Generalized the value range for vector functions. --- lib/data_state.cpp | 265 ++++++++++++++++++++++++++++++++++++++++++- lib/data_state.hpp | 26 ++++- lib/vssolution.cpp | 23 +++- lib/vssolution3d.cpp | 21 ++-- lib/vsvector.cpp | 40 +++++++ lib/vsvector.hpp | 3 + lib/window.cpp | 18 --- 7 files changed, 356 insertions(+), 40 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index d30c1ea9..5f830ed1 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -639,16 +639,12 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) r = rph, i = iph; }; - cmplx_mag_max = 0.; - if (cgrid_f->FESpace()->FEColl()->GetMapType(dim) == FiniteElement::VALUE) { for (int i = 0; i < gf->Size(); i++) { double rval = cgrid_f->real()(i); double ival = cgrid_f->imag()(i); - double mag = hypot(rval, ival); - cmplx_mag_max = max(cmplx_mag_max, mag); rot_ph(rval, ival); (*gf)(i) = cmplx_2_scalar(rval, ival); } @@ -708,8 +704,6 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) { double rval = r_vals(n,d) * w; double ival = i_vals(n,d) * w; - double mag = hypot(rval, ival); - cmplx_mag_max = max(cmplx_mag_max, mag); rot_ph(rval, ival); const double zval = cmplx_2_scalar(rval, ival); z_vals(n,d) = (w != 0.)?(zval / w):(0.); @@ -919,6 +913,265 @@ void DataState::ProjectVectorFEGridFunction() } } +void DataState::FindComplexValueRange(double &minv, double &maxv, + std::function 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 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 vec2scalar, + std::function 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; + } + + // detect common symmetries + + enum class Symmetry { Unknown, Norm, Phase, Linear }; + + Symmetry sym = Symmetry::Unknown; + + { + double a, b; + a = vec2scalar(-1., -1.); + b = vec2scalar(+1., +1.); + 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) + { + minv = (sym == Symmetry::Norm)?(0.):(vec2scalar(-M_PI, -M_PI)); + maxv = vec2scalar(+M_PI, +M_PI); + return; + } + + // real or imag components + + function kernel; + if (sym == Symmetry::Norm) + { + kernel = [&](double r_x, double i_x, double r_y, double i_y) + { + double m_x = hypot(r_x, i_x); + double m_y = hypot(r_y, i_y); + double scal = vec2scalar(m_x, m_y); + if (scale) { scal = scale(scal); } + minv = min(minv, scal); + maxv = max(maxv, scal); + }; + } + else + { + kernel = [&](double r_x, double i_x, double r_y, double i_y) + { + double m_x = hypot(r_x, i_x); + double m_y = hypot(r_y, i_y); + double scal = vec2scalar(m_x, m_y); + if (scale) { scal = scale(scal); } + maxv = max(maxv, scal); + }; + } + + minv = INFINITY; + maxv = 0.; + + const FiniteElementSpace *fes = cgrid_f->FESpace(); + const Mesh *msh = fes->GetMesh(); + const int dim = msh->Dimension(); + MFEM_VERIFY(fes->GetVectorDim() == 2, "Not a 2D vector function!"); + + if (fes->FEColl()->GetMapType(dim) == FiniteElement::VALUE + && fes->FEColl()->GetRangeType(dim) == FiniteElement::SCALAR) + { + const int ndofs = fes->GetNDofs(); + for (int i = 0; i < ndofs; i++) + { + const int ix = fes->DofToVDof(i, 0); + const int iy = fes->DofToVDof(i, 1); + double rval_x = cgrid_f->real()(ix); + double ival_x = cgrid_f->imag()(ix); + double rval_y = cgrid_f->real()(iy); + double ival_y = cgrid_f->imag()(iy); + kernel(rval_x, ival_x, rval_y, ival_y); + } + } + else + { + Array vdofs; + ElementTransformation *Tr; + Vector r_data, i_data; + DenseMatrix r_vals, i_vals; + Vector r_vec, i_vec; + DenseMatrix vshape; + Vector shape; + + const int vdim = fes->GetVDim(); + const int sdim = msh->SpaceDimension(); + + 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(); + + if (fe->GetRangeType() == FiniteElement::SCALAR) + { + r_vals.Reset(r_data.GetData(), nnp, vdim); + i_vals.Reset(i_data.GetData(), nnp, vdim); + 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_x = r_vals(n,0) * w; + double ival_x = i_vals(n,0) * w; + double rval_y = r_vals(n,1) * w; + double ival_y = i_vals(n,1) * w; + kernel(rval_x, ival_x, rval_y, ival_y); + } + } + else + { + vshape.SetSize(nnp, sdim); + for (int n = 0; n < nnp; n++) + { + const IntegrationPoint &ip = ir.IntPoint(n); + Tr->SetIntPoint(&ip); + fe->CalcPhysVShape(*Tr, vshape); + double rval_x = r_data(n) * vshape(n,0); + double ival_x = i_data(n) * vshape(n,0); + double rval_y = r_data(n) * vshape(n,1); + double ival_y = i_data(n) * vshape(n,1); + kernel(rval_x, ival_x, rval_y, ival_y); + } + } + } + } + + if (sym != Symmetry::Norm) { minv = -maxv; } +} + void DataState::ComputeDofsOffsets(std::vector &fespaces) { diff --git a/lib/data_state.hpp b/lib/data_state.hpp index a93ab637..88296cba 100644 --- a/lib/data_state.hpp +++ b/lib/data_state.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include @@ -109,6 +110,13 @@ struct DataState static std::unique_ptr ProjectVectorFEGridFunction(std::unique_ptr gf); + void FindComplexValueRange(double &minv, double &maxv, + std::function = {}) const; + + void FindComplexValueRange(double &minv, double &maxv, + std::function, + std::function = {}) const; + /// Compute the dofs offsets from the grid function vector void ComputeDofsOffsets(std::vector &fespaces); @@ -130,8 +138,7 @@ struct DataState bool fix_elem_orient{false}; bool save_coloring{false}; bool keep_attr{false}; - double cmplx_phase{}; - double cmplx_mag_max{0.}; + double cmplx_phase{0.}; DataState() = default; DataState(DataState &&ss) { *this = std::move(ss); } @@ -257,6 +264,21 @@ struct DataState // or Raviart-Thomas space) with a discontinuous piece-wise polynomial Cartesian // product vector grid function of the same order. 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 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 vec2scal, + std::function scale = {}) const + { if (cgrid_f) { FindComplexValueRange(minv, maxv, vec2scal, scale); } else { minv = maxv = 0.; } } }; #endif // GLVIS_DATA_STATE_HPP diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index f8482c93..0889c73b 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -1153,8 +1153,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(); @@ -1167,11 +1173,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 42a0d7cf..9b82b03f 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -1065,15 +1065,20 @@ void VisualizationSceneSolution3d::FindNewValueRange(bool prepare) GridF->FESpace()->FEColl()->GetMapType(mesh->Dimension()) : FiniteElement::VALUE; - if (shading < Shading::Noncomforming || map_type != (int)FiniteElement::VALUE) - { - minv = sol->Min(); - maxv = sol->Max(); - } - else + win.data_state.FindValueRange(minv, maxv); + + if (minv == 0. && maxv == 0.) { - minv = GridF->Min(); - maxv = GridF->Max(); + 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 39937613..3e201d8d 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -912,6 +912,46 @@ 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, Vec2Scalar, + [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, Vec2Scalar, + [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/window.cpp b/lib/window.cpp index b9a79ef0..bee191e8 100644 --- a/lib/window.cpp +++ b/lib/window.cpp @@ -197,25 +197,7 @@ void Window::GLVisStartVis() void Window::SwitchComplexSolution(DataState::ComplexSolution cmplx_type) { data_state.SwitchComplexSolution(cmplx_type); - auto as = (cmplx_type == DataState::ComplexSolution::Magnitude ? - VisualizationSceneScalarData::Autoscale::MeshAndValue : - VisualizationSceneScalarData::Autoscale::Mesh); - vs->SetAutoscale(as, false); ResetMeshAndSolution(data_state); - switch (cmplx_type) - { - case DataState::ComplexSolution::Magnitude: - break; - case DataState::ComplexSolution::Phase: - vs->SetValueRange(-M_PI, +M_PI); - break; - case DataState::ComplexSolution::Real: - case DataState::ComplexSolution::Imag: - vs->SetValueRange(-data_state.cmplx_mag_max, +data_state.cmplx_mag_max); - break; - default: - std::cerr << "Unknown representation" << std::endl; - } } void Window::SwitchQuadSolution(DataState::QuadSolution quad_type) From 57af8e571344eaeca20540605c63c9cdc6d10d47 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 27 Jan 2026 13:36:57 -0800 Subject: [PATCH 33/44] Refactored scalar function representation in 3D. --- lib/vsvector3d.cpp | 23 +++++++++++++---------- lib/vsvector3d.hpp | 14 +++++++++++++- 2 files changed, 26 insertions(+), 11 deletions(-) diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index 5829561f..d801808d 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -277,8 +277,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() { @@ -286,7 +286,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) + @@ -307,7 +307,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++) @@ -315,7 +315,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++) @@ -323,7 +323,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++) @@ -331,14 +331,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); } @@ -387,7 +390,7 @@ void VisualizationSceneVector3d::Init() drawdisp = 0; drawvector = 0; - scal_func = 0; + scal_func = ScalarFunction::Magnitude; ianim = ianimd = 0; ianimmax = 10; diff --git a/lib/vsvector3d.hpp b/lib/vsvector3d.hpp index 9e2e83a2..35d3cfd4 100644 --- a/lib/vsvector3d.hpp +++ b/lib/vsvector3d.hpp @@ -21,7 +21,19 @@ class VisualizationSceneVector3d : public VisualizationSceneSolution3d protected: mfem::Vector *solx, *soly, *solz; - int drawvector, scal_func; + int drawvector; + enum class ScalarFunction + { + 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; From ac93a9d3f92a14898ab8ebb3594d55feaf497e36 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 27 Jan 2026 14:49:19 -0800 Subject: [PATCH 34/44] Generalized the value range detection to 3D. --- lib/data_state.cpp | 89 ++++++++++++++++++++++---------------------- lib/data_state.hpp | 4 +- lib/vssolution3d.cpp | 8 ++-- lib/vsvector.cpp | 6 ++- lib/vsvector3d.cpp | 46 +++++++++++++++++++++++ lib/vsvector3d.hpp | 2 + 6 files changed, 103 insertions(+), 52 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index 5f830ed1..b3f08b34 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -657,7 +657,6 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) ElementTransformation *Tr; Vector r_data, i_data, z_data; DenseMatrix r_vals, i_vals, z_vals; - Vector r_vec, i_vec, z_vec; DenseMatrix vshape; Vector shape; for (int z = 0; z < msh->GetNE(); z++) @@ -1001,7 +1000,7 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, } void DataState::FindComplexValueRange(double &minv, double &maxv, - std::function vec2scalar, + std::function vec2scalar, std::function scale) const { if (!cgrid_f) @@ -1017,6 +1016,8 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, return; } + const int dim = mesh->Dimension(); + // detect common symmetries enum class Symmetry { Unknown, Norm, Phase, Linear }; @@ -1025,8 +1026,11 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, { double a, b; - a = vec2scalar(-1., -1.); - b = vec2scalar(+1., +1.); + Vector va(dim), vb(dim); + va = -1.; + vb = +1.; + a = vec2scalar(va); + b = vec2scalar(vb); if (a != 0. && b != 0.) { if (a == b) { sym = Symmetry::Norm; } @@ -1053,21 +1057,24 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, // complex phase if (cmplx_sol == ComplexSolution::Phase) { - minv = (sym == Symmetry::Norm)?(0.):(vec2scalar(-M_PI, -M_PI)); - maxv = vec2scalar(+M_PI, +M_PI); + Vector va(dim), vb(dim); + va = -M_PI; + vb = +M_PI; + minv = (sym == Symmetry::Norm)?(0.):(vec2scalar(va)); + maxv = vec2scalar(vb); return; } // real or imag components - function kernel; + function kernel; if (sym == Symmetry::Norm) { - kernel = [&](double r_x, double i_x, double r_y, double i_y) + kernel = [&](const Vector &rd, const Vector &id) { - double m_x = hypot(r_x, i_x); - double m_y = hypot(r_y, i_y); - double scal = vec2scalar(m_x, m_y); + Vector md(dim); + for (int d = 0; d < dim; 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); @@ -1075,11 +1082,11 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, } else { - kernel = [&](double r_x, double i_x, double r_y, double i_y) + kernel = [&](const Vector &rd, const Vector &id) { - double m_x = hypot(r_x, i_x); - double m_y = hypot(r_y, i_y); - double scal = vec2scalar(m_x, m_y); + Vector md(dim); + for (int d = 0; d < dim; d++) { md(d) = hypot(rd(d), id(d)); } + double scal = vec2scalar(md); if (scale) { scal = scale(scal); } maxv = max(maxv, scal); }; @@ -1089,23 +1096,21 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, maxv = 0.; const FiniteElementSpace *fes = cgrid_f->FESpace(); - const Mesh *msh = fes->GetMesh(); - const int dim = msh->Dimension(); - MFEM_VERIFY(fes->GetVectorDim() == 2, "Not a 2D vector function!"); if (fes->FEColl()->GetMapType(dim) == FiniteElement::VALUE && fes->FEColl()->GetRangeType(dim) == FiniteElement::SCALAR) { const int ndofs = fes->GetNDofs(); + Array vdofs; + Vector rd, id; for (int i = 0; i < ndofs; i++) { - const int ix = fes->DofToVDof(i, 0); - const int iy = fes->DofToVDof(i, 1); - double rval_x = cgrid_f->real()(ix); - double ival_x = cgrid_f->imag()(ix); - double rval_y = cgrid_f->real()(iy); - double ival_y = cgrid_f->imag()(iy); - kernel(rval_x, ival_x, rval_y, ival_y); + 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 @@ -1113,15 +1118,11 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, Array vdofs; ElementTransformation *Tr; Vector r_data, i_data; - DenseMatrix r_vals, i_vals; - Vector r_vec, i_vec; + Vector r_vec(dim), i_vec(dim); DenseMatrix vshape; Vector shape; - const int vdim = fes->GetVDim(); - const int sdim = msh->SpaceDimension(); - - for (int z = 0; z < msh->GetNE(); z++) + for (int z = 0; z < mesh->GetNE(); z++) { fes->GetElementVDofs(z, vdofs); cgrid_f->real().GetSubVector(vdofs, r_data); @@ -1134,8 +1135,6 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, if (fe->GetRangeType() == FiniteElement::SCALAR) { - r_vals.Reset(r_data.GetData(), nnp, vdim); - i_vals.Reset(i_data.GetData(), nnp, vdim); shape.SetSize(nnp); for (int n = 0; n < nnp; n++) @@ -1144,26 +1143,28 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, Tr->SetIntPoint(&ip); fe->CalcPhysShape(*Tr, shape); double w = shape(n); - double rval_x = r_vals(n,0) * w; - double ival_x = i_vals(n,0) * w; - double rval_y = r_vals(n,1) * w; - double ival_y = i_vals(n,1) * w; - kernel(rval_x, ival_x, rval_y, ival_y); + for (int d = 0; d < dim; 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, sdim); + vshape.SetSize(nnp, dim); for (int n = 0; n < nnp; n++) { const IntegrationPoint &ip = ir.IntPoint(n); Tr->SetIntPoint(&ip); fe->CalcPhysVShape(*Tr, vshape); - double rval_x = r_data(n) * vshape(n,0); - double ival_x = i_data(n) * vshape(n,0); - double rval_y = r_data(n) * vshape(n,1); - double ival_y = i_data(n) * vshape(n,1); - kernel(rval_x, ival_x, rval_y, ival_y); + for (int d = 0; d < dim; d++) + { + r_vec(d) = r_data(n) * vshape(n,d); + i_vec(d) = i_data(n) * vshape(n,d); + } + kernel(r_vec, i_vec); } } } diff --git a/lib/data_state.hpp b/lib/data_state.hpp index 88296cba..7bf02c59 100644 --- a/lib/data_state.hpp +++ b/lib/data_state.hpp @@ -114,7 +114,7 @@ struct DataState std::function = {}) const; void FindComplexValueRange(double &minv, double &maxv, - std::function, + std::function, std::function = {}) const; /// Compute the dofs offsets from the grid function vector @@ -276,7 +276,7 @@ struct DataState /** 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 vec2scal, + std::function vec2scal, std::function scale = {}) const { if (cgrid_f) { FindComplexValueRange(minv, maxv, vec2scal, scale); } else { minv = maxv = 0.; } } }; diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 9b82b03f..9ad64a48 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -1061,14 +1061,14 @@ 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 (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(); diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index 3e201d8d..7b3d69f7 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -916,7 +916,8 @@ void VisualizationSceneVector::FindNewBox(bool prepare) { VisualizationSceneSolution::FindNewBox(bb.x, bb.y, bb.z); - win.data_state.FindValueRange(minv, maxv, Vec2Scalar, + 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.) @@ -936,7 +937,8 @@ void VisualizationSceneVector::FindNewBox(bool prepare) void VisualizationSceneVector::FindNewValueRange(bool prepare) { - win.data_state.FindValueRange(minv, maxv, Vec2Scalar, + 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.) diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index d801808d..d2cf3699 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -560,6 +560,52 @@ void VisualizationSceneVector3d::NewMeshAndSolution( PrepareDisplacedMesh(); } +void VisualizationSceneVector3d::FindNewValueRange(bool prepare) +{ + function 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 35d3cfd4..0278717c 100644 --- a/lib/vsvector3d.hpp +++ b/lib/vsvector3d.hpp @@ -59,6 +59,8 @@ class VisualizationSceneVector3d : public VisualizationSceneSolution3d void NewMeshAndSolution(const DataState &s) override; + void FindNewValueRange(bool prepare) override; + virtual ~VisualizationSceneVector3d(); std::string GetHelpString() const override; From 660ea8fab2c6a31bbbfc35e38451eb44e52eee85 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 27 Jan 2026 15:39:13 -0800 Subject: [PATCH 35/44] Minor consistency. --- lib/vsvector3d.cpp | 4 ++-- lib/vsvector3d.hpp | 7 +++++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index d2cf3699..24c7d576 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -277,7 +277,7 @@ void VisualizationSceneVector3d::ToggleVectorField(int i) PrepareVectorField(); } -const char *VisualizationSceneVector3d::scal_func_name[(int)ScalarFunction::MAX] +const char *VisualizationSceneVector3d::scal_func_name[(int)ScalarFunction::Max] = {"magnitude", "x-component", "y-component", "z-component"}; void VisualizationSceneVector3d::SetScalarFunction() @@ -340,7 +340,7 @@ void VisualizationSceneVector3d::SetScalarFunction() void VisualizationSceneVector3d::ToggleScalarFunction() { - scal_func = (ScalarFunction)(((int)scal_func + 1) % ((int)ScalarFunction::MAX)); + scal_func = (ScalarFunction)(((int)scal_func + 1) % ((int)ScalarFunction::Max)); cout << "Displaying " << scal_func_name[(int)scal_func] << endl; SetScalarFunction(); FindNewValueRange(true); diff --git a/lib/vsvector3d.hpp b/lib/vsvector3d.hpp index 0278717c..346f117b 100644 --- a/lib/vsvector3d.hpp +++ b/lib/vsvector3d.hpp @@ -22,18 +22,21 @@ class VisualizationSceneVector3d : public VisualizationSceneSolution3d mfem::Vector *solx, *soly, *solz; int drawvector; + enum class ScalarFunction { - MIN = -1, + Invalid = -1, + Min = -1, //---------- Magnitude, Component_X, Component_Y, Component_Z, //---------- - MAX + Max } scal_func; static const char *scal_func_name[]; + double mesh_volume; gl3::GlDrawable vector_buf; gl3::GlDrawable displine_buf; From c3495b4f380d73c04375edbbbced9f50572b2bd5 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 27 Jan 2026 16:10:55 -0800 Subject: [PATCH 36/44] Changed the help string to conditional for quads/complex. --- lib/vssolution.cpp | 25 +++++++++++++++++-------- lib/vssolution3d.cpp | 25 +++++++++++++++++-------- lib/vsvector.cpp | 25 +++++++++++++++++-------- lib/vsvector3d.cpp | 25 +++++++++++++++++-------- 4 files changed, 68 insertions(+), 32 deletions(-) diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 0889c73b..dfbe6341 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -65,10 +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 - << "| Q - Cycle complex 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 @@ -100,10 +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 - << "| Alt+. - Start/stop phase anim. |" << endl - << "| Alt+0/Enter - Phase anim. speed |" << 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 diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 9ad64a48..8e50ab73 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -64,10 +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 - << "| Q - Cycle complex 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 @@ -103,10 +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 - << "| Alt+. - Start/stop phase anim. |" << endl - << "| Alt+0/Enter - Phase anim. speed |" << 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 diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index 7b3d69f7..20cc1e31 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -50,10 +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 - << "| Q - Cycle complex 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 @@ -85,10 +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 - << "| Alt+. - Start/stop phase anim. |" << endl - << "| Alt+0/Enter - Phase anim. speed |" << 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 diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index 24c7d576..f3014567 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -59,10 +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 - << "| Q - Cycle complex 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 @@ -96,10 +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 - << "| Alt+. - Start/stop phase anim. |" << endl - << "| Alt+0/Enter - Phase anim. speed |" << 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 From 85b22b273a34feeebaf8b0eefdfc3a2ce55c22ca Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 28 Jan 2026 11:05:14 -0800 Subject: [PATCH 37/44] Change magnitude/phase to L2 projection. --- lib/data_state.cpp | 167 ++++++++++++++++++++++++++++++++------------- 1 file changed, 118 insertions(+), 49 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index b3f08b34..de6c93a4 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -595,24 +595,29 @@ 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; @@ -625,10 +630,13 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) if (cmplx_type != ComplexSolution::Magnitude) { cout << " (+ animated harmonic phase)"; } cout << endl; } - GridFunction *gf = new GridFunction(cgrid_f->FESpace()); - const FiniteElementSpace *fes = cgrid_f->FESpace(); - const Mesh *msh = fes->GetMesh(); - const int dim = msh->Dimension(); + + 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); @@ -639,9 +647,23 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) r = rph, i = iph; }; - if (cgrid_f->FESpace()->FEColl()->GetMapType(dim) == FiniteElement::VALUE) + GridFunction *gf{}; + if (scal_fes || linear) { - for (int i = 0; i < gf->Size(); i++) + 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); @@ -651,65 +673,112 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) } else { - const int vdim = fes->GetVDim(); - const int sdim = fes->GetVectorDim(); Array vdofs; - ElementTransformation *Tr; Vector r_data, i_data, z_data; DenseMatrix r_vals, i_vals, z_vals; - DenseMatrix vshape; - 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); - z_data.SetSize(vdofs.Size()); - Tr = fes->GetElementTransformation(z); - const FiniteElement *fe = fes->GetFE(z); - const IntegrationRule &ir = fe->GetNodes(); - const int nnp = ir.GetNPoints(); - r_vals.Reset(r_data.GetData(), nnp, vdim); - i_vals.Reset(i_data.GetData(), nnp, vdim); - z_vals.Reset(z_data.GetData(), nnp, vdim); + if (linear) + { + // pointwise projection at DOFs + DenseMatrix vshape; + Vector shape; - if (fe->GetRangeType() == FiniteElement::SCALAR) - { - shape.SetSize(nnp); - } - else + auto kernel = [&](double rdof, double idof, double w) { - vshape.SetSize(nnp, sdim); - } - for (int n = 0; n < nnp; n++) + 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++) { - const IntegrationPoint &ip = ir.IntPoint(n); - Tr->SetIntPoint(&ip); - double w; + 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) { - fe->CalcPhysShape(*Tr, shape); - w = shape(n); + 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 { - fe->CalcPhysVShape(*Tr, vshape); - Vector vec(sdim); - vshape.GetRow(n, vec); - w = vec.Norml2(); + 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); + } } - for (int d = 0; d < vdim; d++) + + 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); + cgrid_f->imag().GetSubVector(vdofs, i_data); + + 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++) { - double rval = r_vals(n,d) * w; - double ival = i_vals(n,d) * w; - rot_ph(rval, ival); - const double zval = cmplx_2_scalar(rval, ival); - z_vals(n,d) = (w != 0.)?(zval / w):(0.); + 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); + gf->SetSubVector(vdofs, z_data); + } } } internal.grid_f.reset(gf); From 4d5410f66aa6f217faeda03c0eb671f2816672cd Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 28 Jan 2026 15:59:29 -0800 Subject: [PATCH 38/44] Minor fixes. --- lib/data_state.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index de6c93a4..cbb7202e 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -752,7 +752,6 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) for (int z = 0; z < mesh->GetNE(); z++) { l2_fes->GetElementVDofs(z, vdofs); - cgrid_f->imag().GetSubVector(vdofs, i_data); ElementTransformation *Tr = fes->GetElementTransformation(z); const FiniteElement *fe = l2_fes->GetFE(z); @@ -1086,6 +1085,7 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, } const int dim = mesh->Dimension(); + const int vdim = cgrid_f->VectorDim(); // detect common symmetries @@ -1095,7 +1095,7 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, { double a, b; - Vector va(dim), vb(dim); + Vector va(vdim), vb(vdim); va = -1.; vb = +1.; a = vec2scalar(va); @@ -1126,7 +1126,7 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, // complex phase if (cmplx_sol == ComplexSolution::Phase) { - Vector va(dim), vb(dim); + Vector va(vdim), vb(vdim); va = -M_PI; vb = +M_PI; minv = (sym == Symmetry::Norm)?(0.):(vec2scalar(va)); @@ -1141,8 +1141,8 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, { kernel = [&](const Vector &rd, const Vector &id) { - Vector md(dim); - for (int d = 0; d < dim; d++) { md(d) = hypot(rd(d), id(d)); } + 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); @@ -1153,8 +1153,8 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, { kernel = [&](const Vector &rd, const Vector &id) { - Vector md(dim); - for (int d = 0; d < dim; d++) { md(d) = hypot(rd(d), id(d)); } + 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); @@ -1187,7 +1187,7 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, Array vdofs; ElementTransformation *Tr; Vector r_data, i_data; - Vector r_vec(dim), i_vec(dim); + Vector r_vec(vdim), i_vec(vdim); DenseMatrix vshape; Vector shape; @@ -1212,7 +1212,7 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, Tr->SetIntPoint(&ip); fe->CalcPhysShape(*Tr, shape); double w = shape(n); - for (int d = 0; d < dim; d++) + 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; @@ -1222,13 +1222,13 @@ void DataState::FindComplexValueRange(double &minv, double &maxv, } else { - vshape.SetSize(nnp, dim); + 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 < dim; d++) + 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); From eec68a74ee5dc29c7c9d7d4c58c19282f422bb1f Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 28 Jan 2026 16:01:08 -0800 Subject: [PATCH 39/44] Added extrusion of complex gfs similar to real to handle 2D3V elements. --- lib/data_state.cpp | 59 ++++++++++++++++++++++++++++++++++------------ lib/data_state.hpp | 3 --- lib/window.cpp | 4 ++-- 3 files changed, 46 insertions(+), 20 deletions(-) diff --git a/lib/data_state.cpp b/lib/data_state.cpp index cbb7202e..36d2fedd 100644 --- a/lib/data_state.cpp +++ b/lib/data_state.cpp @@ -350,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) { @@ -394,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 @@ -414,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(fes2d->FEColl())); + if (cgrid_f) + { + ComplexGridFunction *cgf3d = new ComplexGridFunction(fes3d); + (Vector&)*cgf3d = (const Vector&)*cgrid_f; + cgrid_f->MakeOwner(NULL); + cgf3d->MakeOwner(const_cast(fes2d->FEColl())); + SetCmplxGridFunction(cgf3d); + } + else + { + GridFunction *gf3d = new GridFunction(fes3d); + *gf3d = *grid_f; + grid_f->MakeOwner(NULL); + gf3d->MakeOwner(const_cast(fes2d->FEColl())); + SetGridFunction(gf3d); + } - SetGridFunction(gf3d); delete fes2d; SetMesh(mesh3d); } @@ -785,12 +820,6 @@ void DataState::SetComplexSolution(ComplexSolution cmplx_type, bool print) cmplx_sol = cmplx_type; } -void DataState::SwitchComplexSolution(ComplexSolution cmplx_type, bool print) -{ - SetComplexSolution(cmplx_type, print); - ExtrudeMeshAndSolution(); -} - void DataState::SetQuadSolution(QuadSolution quad_type) { // assume identical order @@ -972,7 +1001,7 @@ void DataState::ProjectVectorFEGridFunction() if (cgrid_f) { internal.cgrid_f = ProjectVectorFEGridFunction(std::move(internal.cgrid_f)); - SwitchComplexSolution(GetComplexSolution(), false); + SetComplexSolution(GetComplexSolution(), false); } else { diff --git a/lib/data_state.hpp b/lib/data_state.hpp index 7bf02c59..47a0d5bc 100644 --- a/lib/data_state.hpp +++ b/lib/data_state.hpp @@ -245,9 +245,6 @@ struct DataState void SetComplexSolution(ComplexSolution type = ComplexSolution::Magnitude, bool print = true); - /// Switch the complex function representation - void SwitchComplexSolution(ComplexSolution type, bool print = true); - /// Get the current representation of complex solution inline ComplexSolution GetComplexSolution() const { return cmplx_sol; } diff --git a/lib/window.cpp b/lib/window.cpp index bee191e8..3b850b15 100644 --- a/lib/window.cpp +++ b/lib/window.cpp @@ -196,7 +196,7 @@ void Window::GLVisStartVis() void Window::SwitchComplexSolution(DataState::ComplexSolution cmplx_type) { - data_state.SwitchComplexSolution(cmplx_type); + data_state.SetComplexSolution(cmplx_type); ResetMeshAndSolution(data_state); } @@ -213,7 +213,7 @@ void Window::UpdateComplexPhase(double ph) DataState::ComplexSolution cs = data_state.GetComplexSolution(); // check if magnitude is viewed, which remains the same if (cs == DataState::ComplexSolution::Magnitude) { return; } - data_state.SwitchComplexSolution(cs, false); + data_state.SetComplexSolution(cs, false); // do not autoscale for animation auto as = vs->GetAutoscale(); vs->SetAutoscale(VisualizationSceneScalarData::Autoscale::None); From 36e067e48f7915ac7bd16481d27eca6c6009c288 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 28 Jan 2026 17:47:48 -0800 Subject: [PATCH 40/44] Workaround for recognition of complex streams. --- lib/stream_reader.cpp | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 4ed16925..be34ca7f 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -106,16 +106,13 @@ bool StreamReader::SupportsDataType(const string &data_type) bool StreamReader::CheckStreamIsComplex(std::istream &solin, bool parallel) { - string buff; - solin >> ws; - getline(solin, buff); - solin.unget(); - const size_t len = buff.length(); - for (size_t i = 0; i < len; i++) { solin.putback(buff[len-1-i]); } - mfem::filter_dos(buff); const char *header = (parallel)?("ParComplexGridFunction"): ("ComplexGridFunction"); - return (buff == header); + 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( From 8752abbb8827bcbe4522e5794f8092235849c4bc Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 28 Jan 2026 18:47:13 -0800 Subject: [PATCH 41/44] Fixed writing of streams. --- lib/stream_reader.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index be34ca7f..ebab0451 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -443,7 +443,7 @@ void StreamReader::WriteStream(std::ostream &os) } else if (data.cgrid_f) { - os << "csolution\n"; + os << "solution\n"; data.mesh->Print(os); data.cgrid_f->Save(os); } From 7bc5bcea88f2c83b42a607f9a9abd8119ef45279 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 28 Jan 2026 19:40:38 -0800 Subject: [PATCH 42/44] Added example 22 test. --- tests/CMakeLists.txt | 2 ++ tests/data | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) 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..1d4a2550 160000 --- a/tests/data +++ b/tests/data @@ -1 +1 @@ -Subproject commit bd2a11d9eabf7dcf2d83cffd906b313e56ab7ef0 +Subproject commit 1d4a2550f7078aa8735ce84c454f17b8331094a9 From 9ff2bc7e4264adc50e6046abc1317a2b6c6311c4 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 28 Jan 2026 19:55:38 -0800 Subject: [PATCH 43/44] Updated example 22 test. --- tests/data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/data b/tests/data index 1d4a2550..f1d01ff0 160000 --- a/tests/data +++ b/tests/data @@ -1 +1 @@ -Subproject commit 1d4a2550f7078aa8735ce84c454f17b8331094a9 +Subproject commit f1d01ff03a5e1d2149e15adef506a803100b45ec From c0bd62dc964d0e67982bad7eafd620cba1ef5ece Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 30 Jan 2026 13:06:04 -0800 Subject: [PATCH 44/44] Fixed (hopefully) handling of Q on Mac. --- lib/aux_vis.cpp | 11 ++++++++--- lib/aux_vis.hpp | 1 + lib/window.cpp | 22 ++++++++++++---------- lib/window.hpp | 3 +++ 4 files changed, 24 insertions(+), 13 deletions(-) diff --git a/lib/aux_vis.cpp b/lib/aux_vis.cpp index 8c886ade..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); @@ -1356,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" }; diff --git a/lib/aux_vis.hpp b/lib/aux_vis.hpp index 98f48cb0..ce7ee31e 100644 --- a/lib/aux_vis.hpp +++ b/lib/aux_vis.hpp @@ -68,6 +68,7 @@ void TouchPinch(SDL_MultiGestureEvent & e); void KeyCtrlP(); void KeyS(); +void KeyqPressed(); void KeyQPressed(); void ToggleThreads(); void ThreadsPauseFunc(SDL_Keymod); diff --git a/lib/window.cpp b/lib/window.cpp index 3b850b15..fe5fe435 100644 --- a/lib/window.cpp +++ b/lib/window.cpp @@ -73,15 +73,6 @@ bool Window::GLVisInitVis(StreamCollection input_streams) locwin = this; - if (data_state.cgrid_f) - { - wnd->setOnKeyDown('Q', SwitchComplexSolution); - } - else 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) @@ -252,9 +243,20 @@ 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) diff --git a/lib/window.hpp b/lib/window.hpp index 9e81cc06..3391a2aa 100644 --- a/lib/window.hpp +++ b/lib/window.hpp @@ -79,6 +79,9 @@ 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;