From b9e0c53faf7c104466c96ac5b5ad2351d289c85b Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 27 Jun 2024 14:48:11 -0700 Subject: [PATCH 01/57] Replaced literals by enum for input type. --- glvis.cpp | 47 +++++++++++++++++++++++++++++++---------------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index c8c97d8e..8364c48e 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -61,7 +61,16 @@ thread_local string extra_caption; bool secure = socketstream::secure_default; // Global variables -int input = 1; +enum InputOptions +{ + INPUT_SERVER_MODE = 1, + INPUT_MESH = 2, + INPUT_SCALAR_SOL = 4, + INPUT_VECTOR_SOL = 8, + //... + INPUT_PARALLEL = 256, +}; +int input = INPUT_SERVER_MODE; thread_local StreamState stream_state; thread_local VisualizationSceneScalarData *vs = NULL; extern thread_local GLVisCommand* glvis_command; @@ -1258,16 +1267,16 @@ int main (int argc, char *argv[]) // set options if (mesh_file != string_none) { - input |= 2; + input |= INPUT_MESH; } if (sol_file != string_none) { - input |= 4; + input |= INPUT_SCALAR_SOL; } if (vec_sol_file != string_none) { sol_file = vec_sol_file; - input |= 8; + input |= INPUT_VECTOR_SOL; } if (gfunc_file != string_none) { @@ -1276,7 +1285,7 @@ int main (int argc, char *argv[]) } if (np > 0) { - input |= 256; + input |= INPUT_PARALLEL; } if (arg_keys != string_none) { @@ -1341,8 +1350,14 @@ int main (int argc, char *argv[]) } // print help for wrong input - if (!(input == 1 || input == 3 || input == 7 || input == 11 || input == 259 || - (stream_state.is_gf && (input == 3 || input == 259)))) + if (!(input == INPUT_SERVER_MODE + || input == (INPUT_SERVER_MODE | INPUT_MESH) + || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_SCALAR_SOL) + || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_VECTOR_SOL) + || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_PARALLEL) + || (stream_state.is_gf + && (input == (INPUT_SERVER_MODE | INPUT_MESH) + || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_PARALLEL))))) { cout << "Invalid combination of mesh/solution options!\n\n"; PrintSampleUsage(cout); @@ -1360,7 +1375,7 @@ int main (int argc, char *argv[]) #endif // server mode, read the mesh and the solution from a socket - if (input == 1) + if (input == INPUT_SERVER_MODE) { // Run server in new thread std::thread serverThread{GLVisServer, portnum, save_stream, @@ -1374,7 +1389,7 @@ int main (int argc, char *argv[]) } else // input != 1, non-server mode { - if (input & 256) + if (input & INPUT_PARALLEL) { ReadParallel(np, stream_state); } @@ -1383,8 +1398,8 @@ int main (int argc, char *argv[]) ReadSerial(stream_state); } - bool use_vector_soln = (input & 8); - bool use_soln = (input & 4); + bool use_vector_soln = (input & INPUT_VECTOR_SOL); + bool use_soln = (input & INPUT_SCALAR_SOL); int field_type; if (use_vector_soln) { @@ -1433,7 +1448,7 @@ void ReadSerial(StreamState& state) state.mesh.reset(new Mesh(meshin, 1, 0, state.fix_elem_orient)); - if (state.is_gf || (input & 4) || (input & 8)) + if (state.is_gf || (input & INPUT_SCALAR_SOL) || (input & INPUT_VECTOR_SOL)) { // get the solution from file bool freesolin = false; @@ -1458,14 +1473,14 @@ void ReadSerial(StreamState& state) state.grid_f.reset(new GridFunction(state.mesh.get(), *solin)); SetGridFunction(state); } - else if (input & 4) + else if (input & INPUT_SCALAR_SOL) { // get rid of NetGen's info line char buff[128]; solin->getline(buff,128); state.sol.Load(*solin, state.mesh->GetNV()); } - else if (input & 8) + else if (input & INPUT_VECTOR_SOL) { state.solu.Load(*solin, state.mesh->GetNV()); state.solv.Load(*solin, state.mesh->GetNV()); @@ -1512,11 +1527,11 @@ void SetGridFunction(StreamState& state) if (state.grid_f->VectorDim() == 1) { state.grid_f->GetNodalValues(state.sol); - input |= 4; + input |= INPUT_SCALAR_SOL; } else { - input |= 8; + input |= INPUT_VECTOR_SOL; } } From f54568d8d2d1502f47fe92be91cc4d5dced449b1 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 27 Jun 2024 17:06:57 -0700 Subject: [PATCH 02/57] Added loading of quad functions to StreamState. --- glvis.cpp | 174 +++++++++++++++++++++++++++++++++++++++++- lib/stream_reader.hpp | 2 + 2 files changed, 175 insertions(+), 1 deletion(-) diff --git a/glvis.cpp b/glvis.cpp index 8364c48e..7f0d1053 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -47,9 +47,11 @@ const char *mesh_file = string_none; const char *sol_file = string_none; const char *vec_sol_file = string_none; const char *gfunc_file = string_none; +const char *qfunc_file = string_none; const char *arg_keys = string_none; int pad_digits = 6; int gf_component = -1; +int qf_component = -1; int window_x = 0; // not a command line option int window_y = 0; // not a command line option int window_w = 400; @@ -98,12 +100,18 @@ void ReadSerial(StreamState& state); // choose grid function component and set the input flag void SetGridFunction(StreamState& state); +// choose quadrature function component and set the input flag +void SetQuadFunction(StreamState& state); + // read the mesh and the solution from multiple files void ReadParallel(int np, StreamState& state); int ReadParMeshAndGridFunction(int np, const char *mesh_prefix, const char *sol_prefix, StreamState& state); +int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, + const char *sol_prefix, StreamState& state); + // Visualize the data in the global variables mesh, sol/grid_f, etc // 0 - scalar data, 1 - vector data, 2 - mesh only, (-1) - unknown bool GLVisInitVis(int field_type, StreamCollection input_streams) @@ -1183,6 +1191,11 @@ int main (int argc, char *argv[]) args.AddOption(&gf_component, "-gc", "--grid-function-component", "Select a grid function component, [0-) or" " -1 for all."); + args.AddOption(&qfunc_file, "-q", "--quadrature-function", + "Quadrature function file to visualize."); + args.AddOption(&qf_component, "-qc", "--quadrature-function-component", + "Select a quadrature function component, [0-) or" + " -1 for all."); args.AddOption(&sol_file, "-s", "--scalar-solution", "Scalar solution (vertex values) file to visualize."); args.AddOption(&vec_sol_file, "-v", "--vector-solution", @@ -1283,6 +1296,11 @@ int main (int argc, char *argv[]) sol_file = gfunc_file; stream_state.is_gf = 255; } + if (qfunc_file != string_none) + { + sol_file = qfunc_file; + stream_state.is_qf = 255; + } if (np > 0) { input |= INPUT_PARALLEL; @@ -1356,6 +1374,9 @@ int main (int argc, char *argv[]) || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_VECTOR_SOL) || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_PARALLEL) || (stream_state.is_gf + && (input == (INPUT_SERVER_MODE | INPUT_MESH) + || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_PARALLEL))) + || (stream_state.is_qf && (input == (INPUT_SERVER_MODE | INPUT_MESH) || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_PARALLEL))))) { @@ -1448,7 +1469,8 @@ void ReadSerial(StreamState& state) state.mesh.reset(new Mesh(meshin, 1, 0, state.fix_elem_orient)); - if (state.is_gf || (input & INPUT_SCALAR_SOL) || (input & INPUT_VECTOR_SOL)) + if (state.is_gf || state.is_qf || (input & INPUT_SCALAR_SOL) || + (input & INPUT_VECTOR_SOL)) { // get the solution from file bool freesolin = false; @@ -1473,6 +1495,11 @@ void ReadSerial(StreamState& state) state.grid_f.reset(new GridFunction(state.mesh.get(), *solin)); SetGridFunction(state); } + else if (state.is_qf) + { + state.quad_f.reset(new QuadratureFunction(state.mesh.get(), *solin)); + SetQuadFunction(state); + } else if (input & INPUT_SCALAR_SOL) { // get rid of NetGen's info line @@ -1535,6 +1562,35 @@ void SetGridFunction(StreamState& state) } } +void SetQuadFunction(StreamState& state) +{ + const int vdim = state.quad_f->GetVDim(); + if (qf_component != -1) + { + if (qf_component < 0 || qf_component >= vdim) + { + cerr << "Invalid component " << qf_component << '.' << endl; + exit(1); + } + QuadratureSpaceBase *qspace = state.quad_f->GetSpace(); + QuadratureFunction *new_qf = new QuadratureFunction(qspace); + for (int i = 0; i < new_qf->Size(); i++) + { + (*new_qf)(i) = (*state.quad_f)(i * vdim + qf_component); + } + state.quad_f->SetOwnsSpace(false); + new_qf->SetOwnsSpace(true); + state.quad_f.reset(new_qf); + } + if (vdim == 1) + { + input |= INPUT_SCALAR_SOL; + } + else + { + input |= INPUT_VECTOR_SOL; + } +} void ReadParallel(int np, StreamState& state) { @@ -1549,6 +1605,15 @@ void ReadParallel(int np, StreamState& state) SetGridFunction(state); } } + else if (state.is_qf) + { + read_err = ReadParMeshAndQuadFunction(np, mesh_file, sol_file, state); + + if (!read_err) + { + SetQuadFunction(state); + } + } else { read_err = ReadParMeshAndGridFunction(np, mesh_file, NULL, @@ -1657,3 +1722,110 @@ int ReadParMeshAndGridFunction(int np, const char *mesh_prefix, return read_err; } + +int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, + const char *sol_prefix, + StreamState& state) +{ + state.mesh = NULL; + + // are the solutions bundled together with the mesh files? + bool same_file = false; + if (sol_prefix) + { + same_file = !strcmp(sol_prefix, mesh_prefix); + state.grid_f = NULL; + } + + Array mesh_array(np); + Array qf_array(np); + mesh_array = NULL; + qf_array = NULL; + + int read_err = 0; + for (int p = 0; p < np; p++) + { + ostringstream fname; + fname << mesh_prefix << '.' << setfill('0') << setw(pad_digits) << p; + named_ifgzstream meshfile(fname.str().c_str()); + if (!meshfile) + { + cerr << "Could not open mesh file: " << fname.str() << '!' << endl; + read_err = 1; + break; + } + + mesh_array[p] = new Mesh(meshfile, 1, 0, state.fix_elem_orient); + + if (!state.keep_attr) + { + // set element and boundary attributes to be the processor number + 1 + for (int i = 0; i < mesh_array[p]->GetNE(); i++) + { + mesh_array[p]->GetElement(i)->SetAttribute(p+1); + } + for (int i = 0; i < mesh_array[p]->GetNBE(); i++) + { + mesh_array[p]->GetBdrElement(i)->SetAttribute(p+1); + } + } + + // read the solution + if (sol_prefix) + { + if (!same_file) + { + ostringstream sol_fname; + sol_fname << sol_prefix << '.' << setfill('0') << setw(pad_digits) << p; + ifgzstream solfile(sol_fname.str().c_str()); + if (!solfile) + { + cerr << "Could not open solution file " + << sol_fname.str() << '!' << endl; + read_err = 2; + break; + } + + qf_array[p] = new QuadratureFunction(mesh_array[p], solfile); + } + else // mesh and solution in the same file + { + qf_array[p] = new QuadratureFunction(mesh_array[p], meshfile); + } + } + } + + if (!read_err) + { + // create the combined mesh and gf + state.mesh.reset(new Mesh(mesh_array, np)); + if (sol_prefix) + { + //assume the same vdim + const int vdim = qf_array[0]->GetVDim(); + //assume the same quadrature rule + QuadratureSpace *qspace = new QuadratureSpace(*state.mesh, + qf_array[0]->GetIntRule(0)); + state.quad_f.reset(new QuadratureFunction(qspace, vdim)); + state.quad_f->SetOwnsSpace(true); + real_t *g_data = state.quad_f->GetData(); + for (int p = 0; p < np; p++) + { + const real_t *l_data = qf_array[p]->GetData(); + const int l_size = qf_array[p]->Size(); + MFEM_ASSERT(g_data + l_size >= state.quad_f->GetData() + state.quad_f->Size(), + "Local parts do not fit to the global quadrature function!"); + memcpy(g_data, l_data, l_size * sizeof(real_t)); + g_data += l_size; + } + } + } + + for (int p = 0; p < np; p++) + { + delete qf_array[np-1-p]; + delete mesh_array[np-1-p]; + } + + return read_err; +} diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index a2d58d66..cabc9f4e 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -27,7 +27,9 @@ struct StreamState std::string keys; std::unique_ptr mesh; std::unique_ptr grid_f; + std::unique_ptr quad_f; int is_gf{0}; + int is_qf{0}; bool fix_elem_orient{false}; bool save_coloring{false}; bool keep_attr{false}; From 722a9decb15b694c59cf439e9b2e218de7ceb5b4 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 27 Jun 2024 17:49:41 -0700 Subject: [PATCH 03/57] Added loading of quadratures from stream. --- glvis.cpp | 33 +++++++++++--------------- lib/stream_reader.cpp | 54 +++++++++++++++++++++++++++++++++++++------ lib/stream_reader.hpp | 3 +++ 3 files changed, 63 insertions(+), 27 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index 7f0d1053..b53fee1b 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -1133,9 +1133,18 @@ void GLVisServer(int portnum, bool save_stream, bool fix_elem_orient, { new_session.state.ReadStreams(input_streams); ofs.precision(8); - ofs << "solution\n"; - new_session.state.mesh->Print(ofs); - new_session.state.grid_f->Save(ofs); + if (new_session.state.grid_f) + { + ofs << "solution\n"; + new_session.state.mesh->Print(ofs); + new_session.state.grid_f->Save(ofs); + } + if (new_session.state.quad_f) + { + ofs << "quadrature\n"; + new_session.state.mesh->Print(ofs); + new_session.state.quad_f->Save(ofs); + } } ofs.close(); cout << "Data saved in " << tmp_file << endl; @@ -1801,23 +1810,7 @@ int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, state.mesh.reset(new Mesh(mesh_array, np)); if (sol_prefix) { - //assume the same vdim - const int vdim = qf_array[0]->GetVDim(); - //assume the same quadrature rule - QuadratureSpace *qspace = new QuadratureSpace(*state.mesh, - qf_array[0]->GetIntRule(0)); - state.quad_f.reset(new QuadratureFunction(qspace, vdim)); - state.quad_f->SetOwnsSpace(true); - real_t *g_data = state.quad_f->GetData(); - for (int p = 0; p < np; p++) - { - const real_t *l_data = qf_array[p]->GetData(); - const int l_size = qf_array[p]->Size(); - MFEM_ASSERT(g_data + l_size >= state.quad_f->GetData() + state.quad_f->Size(), - "Local parts do not fit to the global quadrature function!"); - memcpy(g_data, l_data, l_size * sizeof(real_t)); - g_data += l_size; - } + state.CollectQuadratures(qf_array, np); } } diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 29f7568e..9833824e 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -62,6 +62,27 @@ void StreamState::Extrude1DMeshAndSolution() mesh.reset(mesh2d); } +void StreamState::CollectQuadratures(QuadratureFunction *qf_array[], + int npieces) +{ + //assume the same vdim + const int vdim = qf_array[0]->GetVDim(); + //assume the same quadrature rule + QuadratureSpace *qspace = new QuadratureSpace(*mesh, + qf_array[0]->GetIntRule(0)); + quad_f.reset(new QuadratureFunction(qspace, vdim)); + quad_f->SetOwnsSpace(true); + real_t *g_data = quad_f->GetData(); + for (int p = 0; p < npieces; p++) + { + const real_t *l_data = qf_array[p]->GetData(); + const int l_size = qf_array[p]->Size(); + MFEM_ASSERT(g_data + l_size >= quad_f->GetData() + quad_f->Size(), + "Local parts do not fit to the global quadrature function!"); + memcpy(g_data, l_data, l_size * sizeof(real_t)); + g_data += l_size; + } +} void StreamState::SetMeshSolution() { @@ -199,6 +220,12 @@ int StreamState::ReadStream(istream &is, const string &data_type) grid_f.reset(new GridFunction(mesh.get(), is)); field_type = (grid_f->VectorDim() == 1) ? 0 : 1; } + else if (data_type == "quadrature") + { + mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); + quad_f.reset(new QuadratureFunction(mesh.get(), is)); + field_type = (quad_f->GetVDim() == 1) ? 0 : 1; + } else if (data_type == "mesh") { mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); @@ -340,10 +367,12 @@ int StreamState::ReadStreams(const StreamCollection& input_streams) const int nproc = input_streams.size(); Array mesh_array(nproc); Array gf_array(nproc); + Array qf_array(nproc); std::string data_type; int gf_count = 0; + int qf_count = 0; int field_type = 0; for (int p = 0; p < nproc; p++) @@ -371,7 +400,12 @@ int StreamState::ReadStreams(const StreamCollection& input_streams) } } gf_array[p] = NULL; - if (data_type != "mesh") + 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++; @@ -381,21 +415,27 @@ int StreamState::ReadStreams(const StreamCollection& input_streams) #endif } - if (gf_count > 0 && gf_count != nproc) + if ((gf_count > 0 && gf_count != nproc) + || (qf_count > 0 && qf_count != nproc)) { mfem_error("Input streams contain a mixture of data types!"); } mesh.reset(new Mesh(mesh_array, nproc)); - if (gf_count == 0) + if (gf_count > 0) { - SetMeshSolution(); - field_type = 2; + grid_f.reset(new GridFunction(mesh.get(), gf_array, nproc)); + field_type = (grid_f->VectorDim() == 1) ? 0 : 1; + } + else if (qf_count > 0) + { + CollectQuadratures(qf_array, nproc); + field_type = (quad_f->GetVDim() == 1) ? 0 : 1; } else { - grid_f.reset(new GridFunction(mesh.get(), gf_array, nproc)); - field_type = (grid_f->VectorDim() == 1) ? 0 : 1; + SetMeshSolution(); + field_type = 2; } for (int p = 0; p < nproc; p++) diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index cabc9f4e..d19e3166 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -37,6 +37,9 @@ struct StreamState /// Helper function for visualizing 1D data void Extrude1DMeshAndSolution(); + /// Helper function to build the quadrature function from pieces + void CollectQuadratures(mfem::QuadratureFunction *qf_array[], int npieces); + /// Set a (checkerboard) solution when only the mesh is given void SetMeshSolution(); From 59ae0ffc698c364e5fcdfd4de6ea1c163db15a4f Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 27 Jun 2024 20:15:43 -0700 Subject: [PATCH 04/57] Simplified input options. --- glvis.cpp | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index b53fee1b..fdde35ee 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -1376,18 +1376,21 @@ int main (int argc, char *argv[]) return 0; } + //turn off the server mode if other options are present + if (input & ~INPUT_SERVER_MODE) { input &= ~INPUT_SERVER_MODE; } + // print help for wrong input if (!(input == INPUT_SERVER_MODE - || input == (INPUT_SERVER_MODE | INPUT_MESH) - || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_SCALAR_SOL) - || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_VECTOR_SOL) - || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_PARALLEL) + || input == (INPUT_MESH) + || input == (INPUT_MESH | INPUT_SCALAR_SOL) + || input == (INPUT_MESH | INPUT_VECTOR_SOL) + || input == (INPUT_MESH | INPUT_PARALLEL) || (stream_state.is_gf - && (input == (INPUT_SERVER_MODE | INPUT_MESH) - || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_PARALLEL))) + && (input == (INPUT_MESH) + || input == (INPUT_MESH | INPUT_PARALLEL))) || (stream_state.is_qf - && (input == (INPUT_SERVER_MODE | INPUT_MESH) - || input == (INPUT_SERVER_MODE | INPUT_MESH | INPUT_PARALLEL))))) + && (input == (INPUT_MESH) + || input == (INPUT_MESH | INPUT_PARALLEL))))) { cout << "Invalid combination of mesh/solution options!\n\n"; PrintSampleUsage(cout); From 70af16f683cdd3e1673bb83684177d91185b368b Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 27 Jun 2024 20:51:34 -0700 Subject: [PATCH 05/57] Added enum class for the field type. --- glvis.cpp | 43 +++++++++++++++++++++++++------------------ lib/stream_reader.cpp | 37 +++++++++++++++++++------------------ lib/stream_reader.hpp | 16 ++++++++++++++-- 3 files changed, 58 insertions(+), 38 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index fdde35ee..373cd9f2 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -81,7 +81,8 @@ thread_local communication_thread *comm_thread = NULL; thread_local GeometryRefiner GLVisGeometryRefiner; const char *window_titles[] = { "GLVis [scalar data]", - "GLVis [vector data]", "GLVis [mesh]" + "GLVis [vector data]", + "GLVis [mesh]" }; istream *script = NULL; int scr_running = 0; @@ -113,16 +114,17 @@ int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, const char *sol_prefix, StreamState& state); // Visualize the data in the global variables mesh, sol/grid_f, etc -// 0 - scalar data, 1 - vector data, 2 - mesh only, (-1) - unknown -bool GLVisInitVis(int field_type, StreamCollection input_streams) +bool GLVisInitVis(StreamState::FieldType field_type, + StreamCollection input_streams) { - if (field_type < 0 || field_type > 2) + if (field_type < StreamState::FieldType::MIN + || field_type > StreamState::FieldType::MAX) { return false; } const char *win_title = (window_title == string_default) ? - window_titles[field_type] : window_title; + window_titles[(int)field_type] : window_title; if (InitVisualization(win_title, window_x, window_y, window_w, window_h)) { @@ -138,7 +140,8 @@ bool GLVisInitVis(int field_type, StreamCollection input_streams) } double mesh_range = -1.0; - if (field_type == 0 || field_type == 2) + if (field_type == StreamState::FieldType::SCALAR + || field_type == StreamState::FieldType::MESH) { if (stream_state.grid_f) { @@ -160,7 +163,7 @@ bool GLVisInitVis(int field_type, StreamCollection input_streams) { vss->SetGridFunction(*stream_state.grid_f); } - if (field_type == 2) + if (field_type == StreamState::FieldType::MESH) { vs->OrthogonalProjection = 1; vs->SetLight(false); @@ -179,7 +182,7 @@ bool GLVisInitVis(int field_type, StreamCollection input_streams) { vss->SetGridFunction(stream_state.grid_f.get()); } - if (field_type == 2) + if (field_type == StreamState::FieldType::MESH) { if (stream_state.mesh->Dimension() == 3) { @@ -197,7 +200,7 @@ bool GLVisInitVis(int field_type, StreamCollection input_streams) vss->ToggleDrawMesh(); } } - if (field_type == 2) + if (field_type == StreamState::FieldType::MESH) { if (stream_state.grid_f) { @@ -209,7 +212,7 @@ bool GLVisInitVis(int field_type, StreamCollection input_streams) } } } - else if (field_type == 1) + else if (field_type == StreamState::FieldType::VECTOR) { if (stream_state.mesh->SpaceDimension() == 2) { @@ -252,7 +255,8 @@ bool GLVisInitVis(int field_type, StreamCollection input_streams) vs->SetValueRange(-mesh_range, mesh_range); vs->SetAutoscale(0); } - if (stream_state.mesh->SpaceDimension() == 2 && field_type == 2) + if (stream_state.mesh->SpaceDimension() == 2 + && field_type == StreamState::FieldType::MESH) { SetVisualizationScene(vs, 2, stream_state.keys.c_str()); } @@ -880,7 +884,9 @@ void PlayScript(istream &scr) { plot_caption = c_plot_caption; } - if (GLVisInitVis((stream_state.grid_f->VectorDim() == 1) ? 0 : 1, {})) + if (GLVisInitVis((stream_state.grid_f->VectorDim() == 1) ? + StreamState::FieldType::SCALAR : StreamState::FieldType::VECTOR, + {})) { GetAppWindow()->setOnKeyDown(SDLK_SPACE, ScriptControl); GLVisStartVis(); @@ -904,7 +910,7 @@ struct Session { StreamCollection input_streams; StreamState state; - int ft = -1; + StreamState::FieldType ft = StreamState::FieldType::UNKNOWN; std::thread handler; Session(bool fix_elem_orient, @@ -914,7 +920,7 @@ struct Session state.save_coloring = save_coloring; } - Session(int other_ft, StreamState other_state) + Session(StreamState::FieldType other_ft, StreamState other_state) : state(std::move(other_state)) , ft(other_ft) { } @@ -927,7 +933,7 @@ struct Session void StartSession() { auto funcThread = - [](StreamState thread_state, int ftype, StreamCollection is) + [](StreamState thread_state, StreamState::FieldType ftype, StreamCollection is) { // Set thread-local stream state stream_state = std::move(thread_state); @@ -1433,14 +1439,15 @@ int main (int argc, char *argv[]) bool use_vector_soln = (input & INPUT_VECTOR_SOL); bool use_soln = (input & INPUT_SCALAR_SOL); - int field_type; + StreamState::FieldType field_type; if (use_vector_soln) { - field_type = 1; + field_type = StreamState::FieldType::VECTOR; } else { - field_type = (use_soln) ? 0 : 2; + field_type = (use_soln) ? StreamState::FieldType::SCALAR + : StreamState::FieldType::MESH; } Session single_session(field_type, std::move(stream_state)); single_session.StartSession(); diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 9833824e..a97d4787 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -138,10 +138,10 @@ void StreamState::SetMeshSolution() // Read the content of an input stream (e.g. from socket/file) -int StreamState::ReadStream(istream &is, const string &data_type) +StreamState::FieldType StreamState::ReadStream(istream &is, + const string &data_type) { - // 0 - scalar data, 1 - vector data, 2 - mesh only, (-1) - unknown - int field_type = 0; + FieldType field_type = FieldType::SCALAR; keys.clear(); if (data_type == "fem2d_data") { @@ -150,7 +150,7 @@ int StreamState::ReadStream(istream &is, const string &data_type) } else if (data_type == "vfem2d_data" || data_type == "vfem2d_data_keys") { - field_type = 1; + field_type = FieldType::VECTOR; mesh.reset(new Mesh(is, 0, 0, fix_elem_orient)); solu.Load(is, mesh->GetNV()); solv.Load(is, mesh->GetNV()); @@ -166,7 +166,7 @@ int StreamState::ReadStream(istream &is, const string &data_type) } else if (data_type == "vfem3d_data" || data_type == "vfem3d_data_keys") { - field_type = 1; + field_type = FieldType::VECTOR; mesh.reset(new Mesh(is, 0, 0, fix_elem_orient)); solu.Load(is, mesh->GetNV()); solv.Load(is, mesh->GetNV()); @@ -187,7 +187,7 @@ int StreamState::ReadStream(istream &is, const string &data_type) } else if (data_type == "vfem2d_gf_data" || data_type == "vfem2d_gf_data_keys") { - field_type = 1; + field_type = FieldType::VECTOR; mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); grid_f.reset(new GridFunction(mesh.get(), is)); if (data_type == "vfem2d_gf_data_keys") @@ -206,7 +206,7 @@ int StreamState::ReadStream(istream &is, const string &data_type) } else if (data_type == "vfem3d_gf_data" || data_type == "vfem3d_gf_data_keys") { - field_type = 1; + field_type = FieldType::VECTOR; mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); grid_f.reset(new GridFunction(mesh.get(), is)); if (data_type == "vfem3d_gf_data_keys") @@ -218,19 +218,19 @@ int StreamState::ReadStream(istream &is, const string &data_type) { mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); grid_f.reset(new GridFunction(mesh.get(), is)); - field_type = (grid_f->VectorDim() == 1) ? 0 : 1; + field_type = (grid_f->VectorDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR; } else if (data_type == "quadrature") { mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); quad_f.reset(new QuadratureFunction(mesh.get(), is)); - field_type = (quad_f->GetVDim() == 1) ? 0 : 1; + field_type = (quad_f->GetVDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR; } else if (data_type == "mesh") { mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); SetMeshSolution(); - field_type = 2; + field_type = FieldType::MESH; } else if (data_type == "raw_scalar_2d") { @@ -345,16 +345,16 @@ int StreamState::ReadStream(istream &is, const string &data_type) delete vertices[i]; } - field_type = 0; + field_type = FieldType::SCALAR; } else { - field_type = -1; + field_type = FieldType::UNKNOWN; cerr << "Unknown data format" << endl; cerr << data_type << endl; } - if (field_type >= 0 && field_type <= 2) + if (field_type > FieldType::MIN && field_type < FieldType::MAX) { Extrude1DMeshAndSolution(); } @@ -362,7 +362,8 @@ int StreamState::ReadStream(istream &is, const string &data_type) return field_type; } -int StreamState::ReadStreams(const StreamCollection& input_streams) +StreamState::FieldType StreamState::ReadStreams(const StreamCollection& + input_streams) { const int nproc = input_streams.size(); Array mesh_array(nproc); @@ -373,7 +374,7 @@ int StreamState::ReadStreams(const StreamCollection& input_streams) int gf_count = 0; int qf_count = 0; - int field_type = 0; + FieldType field_type = FieldType::SCALAR; for (int p = 0; p < nproc; p++) { @@ -425,17 +426,17 @@ int StreamState::ReadStreams(const StreamCollection& input_streams) if (gf_count > 0) { grid_f.reset(new GridFunction(mesh.get(), gf_array, nproc)); - field_type = (grid_f->VectorDim() == 1) ? 0 : 1; + field_type = (grid_f->VectorDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR; } else if (qf_count > 0) { CollectQuadratures(qf_array, nproc); - field_type = (quad_f->GetVDim() == 1) ? 0 : 1; + field_type = (quad_f->GetVDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR; } else { SetMeshSolution(); - field_type = 2; + field_type = FieldType::MESH; } for (int p = 0; p < nproc; p++) diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index d19e3166..f1604724 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -34,6 +34,18 @@ struct StreamState bool save_coloring{false}; bool keep_attr{false}; + enum class FieldType + { + UNKNOWN = -1, + MIN = -1, + //---------- + SCALAR, + VECTOR, + MESH, + //---------- + MAX + }; + /// Helper function for visualizing 1D data void Extrude1DMeshAndSolution(); @@ -43,9 +55,9 @@ struct StreamState /// Set a (checkerboard) solution when only the mesh is given void SetMeshSolution(); - int ReadStream(std::istream &is, const std::string &data_type); + FieldType ReadStream(std::istream &is, const std::string &data_type); - int ReadStreams(const StreamCollection& input_streams); + FieldType ReadStreams(const StreamCollection& input_streams); /// Sets a new mesh and solution from another StreamState object, and /// updates the given VisualizationScene pointer with the new data. From 9f46ae84181e1036462daaf36da1bdd343f31dc9 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 27 Jun 2024 21:46:20 -0700 Subject: [PATCH 06/57] Added extrusion of 1D quad functions to 2D. --- lib/stream_reader.cpp | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index a97d4787..fa070740 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -17,6 +17,10 @@ using namespace std; using namespace mfem; +/// Helper function for extrusion of 1D quadrature functions to 2D +QuadratureFunction* Extrude1DQuadFunction(Mesh *mesh, Mesh *mesh2d, + QuadratureFunction *qf, int ny); + void StreamState::Extrude1DMeshAndSolution() { if (mesh->Dimension() != 1 || mesh->SpaceDimension() != 1) @@ -49,6 +53,13 @@ void StreamState::Extrude1DMeshAndSolution() grid_f.reset(grid_f_2d); } + else if (quad_f) + { + QuadratureFunction *quad_f_2d = + Extrude1DQuadFunction(mesh.get(), mesh2d, quad_f.get(), 1); + + quad_f.reset(quad_f_2d); + } else if (sol.Size() == mesh->GetNV()) { Vector sol2d(mesh2d->GetNV()); @@ -527,3 +538,33 @@ bool StreamState::SetNewMeshAndSolution(StreamState new_state, return false; } } + +QuadratureFunction *Extrude1DQuadFunction(Mesh *mesh, Mesh *mesh2d, + QuadratureFunction *qf, int ny) +{ + //assume identical orders + const int order = qf->GetIntRule(0).GetOrder(); + const int vdim = qf->GetVDim(); + QuadratureSpace *qspace2d = new QuadratureSpace(mesh2d, order); + QuadratureFunction *qf2d = new QuadratureFunction(qspace2d); + qf2d->SetOwnsSpace(true); + + DenseMatrix vals, vals2d; + for (int ix = 0; ix < mesh->GetNE(); ix++) + { + qf->GetValues(ix, vals); + const int nq = vals.Width(); + for (int iy = 0; iy < ny; iy++) + { + qf2d->GetValues(ix*ny+iy, vals2d); + for (int qx = 0; qx < nq; qx++) + for (int qy = 0; qy < nq; qy++) + for (int d = 0; d < vdim; d++) + { + vals2d(d, qy*nq+qx) = vals(d, qx); + } + } + } + + return qf2d; +} From 47a2a1da6a5d3b36fb292d2427959b6dd88700cc Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 27 Jun 2024 22:11:36 -0700 Subject: [PATCH 07/57] Minor reordering to priritize the original quad function. --- glvis.cpp | 12 ++++++------ lib/stream_reader.cpp | 16 ++++++++-------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index 373cd9f2..408f0f03 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -1139,18 +1139,18 @@ void GLVisServer(int portnum, bool save_stream, bool fix_elem_orient, { new_session.state.ReadStreams(input_streams); ofs.precision(8); - if (new_session.state.grid_f) - { - ofs << "solution\n"; - new_session.state.mesh->Print(ofs); - new_session.state.grid_f->Save(ofs); - } if (new_session.state.quad_f) { ofs << "quadrature\n"; new_session.state.mesh->Print(ofs); new_session.state.quad_f->Save(ofs); } + else if (new_session.state.grid_f) + { + ofs << "solution\n"; + new_session.state.mesh->Print(ofs); + new_session.state.grid_f->Save(ofs); + } } ofs.close(); cout << "Data saved in " << tmp_file << endl; diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index fa070740..b95fd4be 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -46,20 +46,20 @@ void StreamState::Extrude1DMeshAndSolution() Mesh *mesh2d = Extrude1D(mesh.get(), 1, 0.1*(xmax - xmin)); - if (grid_f) - { - GridFunction *grid_f_2d = - Extrude1DGridFunction(mesh.get(), mesh2d, grid_f.get(), 1); - - grid_f.reset(grid_f_2d); - } - else if (quad_f) + if (quad_f) { QuadratureFunction *quad_f_2d = Extrude1DQuadFunction(mesh.get(), mesh2d, quad_f.get(), 1); quad_f.reset(quad_f_2d); } + else if (grid_f) + { + GridFunction *grid_f_2d = + Extrude1DGridFunction(mesh.get(), mesh2d, grid_f.get(), 1); + + grid_f.reset(grid_f_2d); + } else if (sol.Size() == mesh->GetNV()) { Vector sol2d(mesh2d->GetNV()); From dd7bd9651bf0ff23dda2ee451a45299f3296be25 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 27 Jun 2024 23:18:53 -0700 Subject: [PATCH 08/57] Added conversion to high-order L2 grid function. --- glvis.cpp | 2 ++ lib/stream_reader.cpp | 17 +++++++++++++++++ lib/stream_reader.hpp | 9 +++++++++ 3 files changed, 28 insertions(+) diff --git a/glvis.cpp b/glvis.cpp index 408f0f03..3695d2e7 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -1609,6 +1609,8 @@ void SetQuadFunction(StreamState& state) { input |= INPUT_VECTOR_SOL; } + + state.SetQuadSolution(); } void ReadParallel(int np, StreamState& state) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index b95fd4be..cf1a7029 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -147,6 +147,21 @@ void StreamState::SetMeshSolution() } } +void StreamState::SetQuadSolution(QuadSolution type) +{ + //if(type == QuadSolution::HO_L2) + { + //assume identical order + const int order = quad_f->GetIntRule(0).GetOrder()/2;//<---Gauss-Legendre + FiniteElementCollection *fec = new L2_FECollection(order, mesh->Dimension()); + FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, + quad_f->GetVDim(), Ordering::byVDIM); + GridFunction *gf = new GridFunction(fes); + gf->MakeOwner(fec); + *gf = *quad_f; + grid_f.reset(gf); + } +} // Read the content of an input stream (e.g. from socket/file) StreamState::FieldType StreamState::ReadStream(istream &is, @@ -235,6 +250,7 @@ StreamState::FieldType StreamState::ReadStream(istream &is, { mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); quad_f.reset(new QuadratureFunction(mesh.get(), is)); + SetQuadSolution(); field_type = (quad_f->GetVDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR; } else if (data_type == "mesh") @@ -442,6 +458,7 @@ StreamState::FieldType StreamState::ReadStreams(const StreamCollection& else if (qf_count > 0) { CollectQuadratures(qf_array, nproc); + SetQuadSolution(); field_type = (quad_f->GetVDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR; } else diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index f1604724..ace50bdc 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -46,6 +46,12 @@ struct StreamState MAX }; + enum class QuadSolution + { + LOR_GLL, + HO_L2, + }; + /// Helper function for visualizing 1D data void Extrude1DMeshAndSolution(); @@ -55,6 +61,9 @@ struct StreamState /// Set a (checkerboard) solution when only the mesh is given void SetMeshSolution(); + /// Set a quadrature function solution producing a proxy grid function + void SetQuadSolution(QuadSolution type = QuadSolution::LOR_GLL); + FieldType ReadStream(std::istream &is, const std::string &data_type); FieldType ReadStreams(const StreamCollection& input_streams); From 8f22e3e0f97ccdc621f14f0a230ca1ac0cf0d072 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 28 Jun 2024 07:13:03 -0700 Subject: [PATCH 09/57] Used reference on quad function. --- lib/stream_reader.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index cf1a7029..3764c796 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -156,9 +156,9 @@ void StreamState::SetQuadSolution(QuadSolution type) FiniteElementCollection *fec = new L2_FECollection(order, mesh->Dimension()); FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, quad_f->GetVDim(), Ordering::byVDIM); - GridFunction *gf = new GridFunction(fes); + MFEM_ASSERT(quad_f->Size() == fes->GetVSize(), "Size mismatch"); + GridFunction *gf = new GridFunction(fes, *quad_f, 0); gf->MakeOwner(fec); - *gf = *quad_f; grid_f.reset(gf); } } From add5839b23a5970631e12f5085185c71a2b3c7a9 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 28 Jun 2024 07:20:31 -0700 Subject: [PATCH 10/57] Fixed extrusion of 1D quad functions. --- lib/stream_reader.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 3764c796..a084046a 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -52,6 +52,7 @@ void StreamState::Extrude1DMeshAndSolution() Extrude1DQuadFunction(mesh.get(), mesh2d, quad_f.get(), 1); quad_f.reset(quad_f_2d); + SetQuadSolution(); } else if (grid_f) { @@ -149,6 +150,12 @@ void StreamState::SetMeshSolution() void StreamState::SetQuadSolution(QuadSolution type) { + /*if(type == QuadSolution::LOR_GLL) + { + //assume identical order + const int order = quad_f->GetIntRule(0).Size(); + Mesh *mesh_lor = new Mesh(Mesh::MakeRefined(*mesh, order, )); + }*/ //if(type == QuadSolution::HO_L2) { //assume identical order From ec2f1134ba4e48d30f08bfbac72bb1cc56b856c4 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 28 Jun 2024 07:22:15 -0700 Subject: [PATCH 11/57] Fixed field type check. --- glvis.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index 3695d2e7..705a9714 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -117,8 +117,8 @@ int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, bool GLVisInitVis(StreamState::FieldType field_type, StreamCollection input_streams) { - if (field_type < StreamState::FieldType::MIN - || field_type > StreamState::FieldType::MAX) + if (field_type <= StreamState::FieldType::MIN + || field_type >= StreamState::FieldType::MAX) { return false; } From c3137e64612d5f20a8590e8a077297573273c121 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 28 Jun 2024 07:33:22 -0700 Subject: [PATCH 12/57] Updated sample usage. --- glvis.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/glvis.cpp b/glvis.cpp index 705a9714..4f03e57a 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -1471,7 +1471,11 @@ void PrintSampleUsage(ostream &os) "Visualize mesh and solution (grid function):\n" " glvis -m -g [-gc ]\n" "Visualize parallel mesh and solution (grid function):\n" - " glvis -np <#proc> -m [-g ]\n\n" + " glvis -np <#proc> -m [-g ]\n" + "Visualize mesh and quadrature function:\n" + " glvis -m -q [-qc ]\n" + "Visualize parallel mesh and quadrature function:\n" + " glvis -np <#proc> -m [-q ]\n\n" "All Options:\n"; } From 3c80af4cb626a7ca674f17c89f8ee3db85981551 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 28 Jun 2024 07:49:34 -0700 Subject: [PATCH 13/57] Added a message about quad func conversion. --- lib/stream_reader.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index a084046a..5b37127f 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -164,6 +164,8 @@ void StreamState::SetQuadSolution(QuadSolution type) FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, quad_f->GetVDim(), Ordering::byVDIM); MFEM_ASSERT(quad_f->Size() == fes->GetVSize(), "Size mismatch"); + cout << "Representing quadrature by grid function of order " + << order << endl; GridFunction *gf = new GridFunction(fes, *quad_f, 0); gf->MakeOwner(fec); grid_f.reset(gf); From 9865d8bac86c35c4abc8b6470d99d4a54bea9577 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 28 Jun 2024 08:31:48 -0700 Subject: [PATCH 14/57] Added LOR visualization of quad functions. Gauss-Lobatto for the moment. --- lib/stream_reader.cpp | 34 ++++++++++++++++++++++++++-------- 1 file changed, 26 insertions(+), 8 deletions(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 5b37127f..34687191 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -150,16 +150,34 @@ void StreamState::SetMeshSolution() void StreamState::SetQuadSolution(QuadSolution type) { - /*if(type == QuadSolution::LOR_GLL) + //assume identical order + const int order = quad_f->GetIntRule(0).GetOrder()/2;//<---Gauss-Legendre + if (type == QuadSolution::LOR_GLL) { - //assume identical order - const int order = quad_f->GetIntRule(0).Size(); - Mesh *mesh_lor = new Mesh(Mesh::MakeRefined(*mesh, order, )); - }*/ - //if(type == QuadSolution::HO_L2) + const int ref_factor = order + 1; + if (ref_factor > 1) + { + Mesh *mesh_lor = new Mesh(Mesh::MakeRefined(*mesh, ref_factor, + BasisType::GaussLobatto)); + FiniteElementCollection *fec = new L2_FECollection(0, mesh->Dimension()); + FiniteElementSpace *fes = new FiniteElementSpace(mesh_lor, fec, + quad_f->GetVDim(), Ordering::byVDIM); + MFEM_ASSERT(quad_f->Size() == fes->GetVSize(), "Size mismatch"); + cout << "Representing quadrature by piecewise-constant function on mesh refined " + << ref_factor << " times" << endl; + GridFunction *gf = new GridFunction(fes, *quad_f, 0); + gf->MakeOwner(fec); + grid_f.reset(gf); + mesh.reset(mesh_lor); + } + else + { + //low-order + type = QuadSolution::HO_L2; + } + } + if (type == QuadSolution::HO_L2) { - //assume identical order - const int order = quad_f->GetIntRule(0).GetOrder()/2;//<---Gauss-Legendre FiniteElementCollection *fec = new L2_FECollection(order, mesh->Dimension()); FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, quad_f->GetVDim(), Ordering::byVDIM); From 02b7a06bedc32238370d6001a464c07e035a5651 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 28 Jun 2024 10:51:00 -0700 Subject: [PATCH 15/57] Moved writing the stream to SteamState. --- glvis.cpp | 14 +------------- lib/stream_reader.cpp | 17 +++++++++++++++++ lib/stream_reader.hpp | 2 ++ 3 files changed, 20 insertions(+), 13 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index 4f03e57a..aec5c998 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -1138,19 +1138,7 @@ void GLVisServer(int portnum, bool save_stream, bool fix_elem_orient, else { new_session.state.ReadStreams(input_streams); - ofs.precision(8); - if (new_session.state.quad_f) - { - ofs << "quadrature\n"; - new_session.state.mesh->Print(ofs); - new_session.state.quad_f->Save(ofs); - } - else if (new_session.state.grid_f) - { - ofs << "solution\n"; - new_session.state.mesh->Print(ofs); - new_session.state.grid_f->Save(ofs); - } + new_session.state.WriteStream(ofs); } ofs.close(); cout << "Data saved in " << tmp_file << endl; diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 34687191..441162d7 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -505,6 +505,23 @@ StreamState::FieldType StreamState::ReadStreams(const StreamCollection& return field_type; } +void StreamState::WriteStream(std::ostream &os) +{ + os.precision(8); + if (quad_f) + { + os << "quadrature\n"; + mesh->Print(os); + quad_f->Save(os); + } + else if (grid_f) + { + os << "solution\n"; + mesh->Print(os); + grid_f->Save(os); + } +} + // 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. diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index ace50bdc..eefce378 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -68,6 +68,8 @@ struct StreamState FieldType ReadStreams(const StreamCollection& input_streams); + void WriteStream(std::ostream &os); + /// Sets a new mesh and solution from another StreamState object, and /// updates the given VisualizationScene pointer with the new data. /// From 23213d540b6aec6e5ef244cea85bc2df28f98883 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 28 Jun 2024 11:30:29 -0700 Subject: [PATCH 16/57] Saved the original mesh before LOR. --- lib/stream_reader.cpp | 16 +++++++++++++++- lib/stream_reader.hpp | 1 + 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 441162d7..defbacf0 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -152,6 +152,12 @@ void StreamState::SetQuadSolution(QuadSolution type) { //assume identical order const int order = quad_f->GetIntRule(0).GetOrder()/2;//<---Gauss-Legendre + //use the original mesh when available + if (mesh_quad.get()) + { + mesh.swap(mesh_quad); + mesh_quad.reset(); + } if (type == QuadSolution::LOR_GLL) { const int ref_factor = order + 1; @@ -168,6 +174,7 @@ void StreamState::SetQuadSolution(QuadSolution type) GridFunction *gf = new GridFunction(fes, *quad_f, 0); gf->MakeOwner(fec); grid_f.reset(gf); + mesh.swap(mesh_quad); mesh.reset(mesh_lor); } else @@ -511,7 +518,14 @@ void StreamState::WriteStream(std::ostream &os) if (quad_f) { os << "quadrature\n"; - mesh->Print(os); + if (mesh_quad.get()) + { + mesh_quad->Print(os); + } + else + { + mesh->Print(os); + } quad_f->Save(os); } else if (grid_f) diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index eefce378..ecaad001 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -26,6 +26,7 @@ struct StreamState mfem::Vector sol, solu, solv, solw, normals; std::string keys; std::unique_ptr mesh; + std::unique_ptr mesh_quad; std::unique_ptr grid_f; std::unique_ptr quad_f; int is_gf{0}; From ae845f63c0ff1090c566fbc73f0f63c4d36196d1 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 28 Jun 2024 11:58:15 -0700 Subject: [PATCH 17/57] Added scripts support for quadratures. --- glvis.cpp | 113 +++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 112 insertions(+), 1 deletion(-) diff --git a/glvis.cpp b/glvis.cpp index aec5c998..de073cf8 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -322,6 +322,46 @@ int ScriptReadSolution(istream &scr, StreamState& state) return 0; } +int ScriptReadQuadrature(istream &scr, StreamState& state) +{ + string mword,sword; + + cout << "Script: quadrature: " << flush; + // read the mesh + scr >> ws >> mword; // mesh filename (can't contain spaces) + cout << "mesh: " << mword << "; " << flush; + named_ifgzstream imesh(mword.c_str()); + if (!imesh) + { + cout << "Can not open mesh file: " << mword << endl; + return 1; + } + state.mesh.reset(new Mesh(imesh, 1, 0, state.fix_elem_orient)); + + // read the quadrature (QuadratureFunction) + scr >> ws >> sword; + if (sword == mword) // mesh and quadrature in the same file + { + cout << "quadrature: " << mword << endl; + state.quad_f.reset(new QuadratureFunction(state.mesh.get(), imesh)); + } + else + { + cout << "quadrature: " << sword << endl; + ifgzstream isol(sword.c_str()); + if (!isol) + { + cout << "Can not open quadrature file: " << sword << endl; + return 2; + } + state.quad_f.reset(new QuadratureFunction(state.mesh.get(), isol)); + } + + state.Extrude1DMeshAndSolution(); + + return 0; +} + int ScriptReadParSolution(istream &scr, StreamState& state) { int np, scr_keep_attr, err_read; @@ -356,6 +396,40 @@ int ScriptReadParSolution(istream &scr, StreamState& state) return err_read; } +int ScriptReadParQuadrature(istream &scr, StreamState& state) +{ + int np, scr_keep_attr, err_read; + string mesh_prefix, quad_prefix; + + cout << "Script: pquadrature: " << flush; + // read number of processors + scr >> np; + cout << "# processors: " << np << "; " << flush; + // read the mesh prefix + scr >> ws >> mesh_prefix; // mesh prefix (can't contain spaces) + cout << "mesh prefix: " << mesh_prefix << "; " << flush; + scr >> ws >> scr_keep_attr; + if (scr_keep_attr) + { + cout << "(real attributes); " << flush; + } + else + { + cout << "(processor attributes); " << flush; + } + // read the quadrature prefix + scr >> ws >> quad_prefix; + cout << "quadrature prefix: " << quad_prefix << endl; + + err_read = ReadParMeshAndQuadFunction(np, mesh_prefix.c_str(), + quad_prefix.c_str(), state); + if (!err_read) + { + state.Extrude1DMeshAndSolution(); + } + return err_read; +} + int ScriptReadDisplMesh(istream &scr, StreamState& state) { StreamState meshstate; @@ -455,7 +529,8 @@ void ExecuteScriptCommand() scr_level = 0; } } - else if (word == "solution" || word == "mesh" || word == "psolution") + else if (word == "solution" || word == "mesh" || word == "psolution" + || word == "quadrature" || word == "pquadrature") { StreamState new_state; @@ -467,6 +542,14 @@ void ExecuteScriptCommand() continue; } } + else if (word == "quadrature") + { + if (ScriptReadQuadrature(scr, new_state)) + { + done_one_command = 1; + continue; + } + } else if (word == "mesh") { if (ScriptReadDisplMesh(scr, new_state)) @@ -489,6 +572,14 @@ void ExecuteScriptCommand() continue; } } + else if (word == "pquadrature") + { + if (ScriptReadParQuadrature(scr, new_state)) + { + done_one_command = 1; + continue; + } + } if (stream_state.SetNewMeshAndSolution(std::move(new_state), vs)) { @@ -843,6 +934,16 @@ void PlayScript(istream &scr) // start the visualization break; } + else if (word == "quadrature") + { + if (ScriptReadQuadrature(scr, stream_state)) + { + return; + } + + // start the visualization + break; + } else if (word == "psolution") { if (ScriptReadParSolution(scr, stream_state)) @@ -853,6 +954,16 @@ void PlayScript(istream &scr) // start the visualization break; } + else if (word == "pquadrature") + { + if (ScriptReadParQuadrature(scr, stream_state)) + { + return; + } + + // start the visualization + break; + } else if (word == "mesh") { if (ScriptReadDisplMesh(scr, stream_state)) From 22f92debd81a719c6a25fbfb7619876e91def5d4 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 28 Jun 2024 12:14:37 -0700 Subject: [PATCH 18/57] Added moving of quadrature and its mesh in StreamState. --- lib/stream_reader.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index defbacf0..5abf7126 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -570,6 +570,8 @@ bool StreamState::SetNewMeshAndSolution(StreamState new_state, { std::unique_ptr new_m = std::move(new_state.mesh); std::unique_ptr new_g = std::move(new_state.grid_f); + std::unique_ptr new_mq = std::move(new_state.mesh_quad); + std::unique_ptr new_q = std::move(new_state.quad_f); if (new_m->SpaceDimension() == 2) { if (new_g->VectorDim() == 1) @@ -606,6 +608,8 @@ bool StreamState::SetNewMeshAndSolution(StreamState new_state, } grid_f = std::move(new_g); mesh = std::move(new_m); + quad_f = std::move(new_q); + mesh_quad = std::move(new_mq); return true; } else From 43b76ed8f041321387c44c31f5ac6479fd5cdd5c Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 1 Jul 2024 10:20:08 -0700 Subject: [PATCH 19/57] Added some resets of the mesh function. --- glvis.cpp | 3 +++ lib/stream_reader.cpp | 2 ++ 2 files changed, 5 insertions(+) diff --git a/glvis.cpp b/glvis.cpp index de073cf8..0747af46 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -337,6 +337,7 @@ int ScriptReadQuadrature(istream &scr, StreamState& state) return 1; } state.mesh.reset(new Mesh(imesh, 1, 0, state.fix_elem_orient)); + state.mesh_quad.reset(); // read the quadrature (QuadratureFunction) scr >> ws >> sword; @@ -1590,6 +1591,7 @@ void ReadSerial(StreamState& state) } state.mesh.reset(new Mesh(meshin, 1, 0, state.fix_elem_orient)); + state.mesh_quad.reset(); if (state.is_gf || state.is_qf || (input & INPUT_SCALAR_SOL) || (input & INPUT_VECTOR_SOL)) @@ -1923,6 +1925,7 @@ int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, { // create the combined mesh and gf state.mesh.reset(new Mesh(mesh_array, np)); + state.mesh_quad.reset(); if (sol_prefix) { state.CollectQuadratures(qf_array, np); diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 5abf7126..a36e7529 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -283,6 +283,7 @@ StreamState::FieldType StreamState::ReadStream(istream &is, else if (data_type == "quadrature") { mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); + mesh_quad.reset(); quad_f.reset(new QuadratureFunction(mesh.get(), is)); SetQuadSolution(); field_type = (quad_f->GetVDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR; @@ -484,6 +485,7 @@ StreamState::FieldType StreamState::ReadStreams(const StreamCollection& } mesh.reset(new Mesh(mesh_array, nproc)); + mesh_quad.reset(); if (gf_count > 0) { grid_f.reset(new GridFunction(mesh.get(), gf_array, nproc)); From 7ae8e3f86dfcbad5b7d9b33543b13469dac7b132 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 1 Jul 2024 10:39:47 -0700 Subject: [PATCH 20/57] Added a command for loading quadrature. --- lib/threads.cpp | 110 ++++++++++++++++++++++++++++++++++++++++++++++-- lib/threads.hpp | 2 + 2 files changed, 109 insertions(+), 3 deletions(-) diff --git a/lib/threads.cpp b/lib/threads.cpp index bf3c128c..4697ce79 100644 --- a/lib/threads.cpp +++ b/lib/threads.cpp @@ -96,6 +96,24 @@ int GLVisCommand::NewMeshAndSolution(std::unique_ptr _new_m, return 0; } +int GLVisCommand::NewMeshAndQuadrature(std::unique_ptr _new_m, + std::unique_ptr _new_q) +{ + if (lock() < 0) + { + return -1; + } + command = NEW_MESH_AND_SOLUTION; + new_state.mesh = std::move(_new_m); + new_state.mesh_quad.reset(); + new_state.quad_f = std::move(_new_q); + if (signal() < 0) + { + return -2; + } + return 0; +} + int GLVisCommand::Screenshot(const char *filename) { if (lock() < 0) @@ -427,9 +445,16 @@ int GLVisCommand::Execute() double mesh_range = -1.0; if (!new_state.grid_f) { - new_state.save_coloring = false; - new_state.SetMeshSolution(); - mesh_range = new_state.grid_f->Max() + 1.0; + if (!new_state.quad_f) + { + new_state.save_coloring = false; + new_state.SetMeshSolution(); + mesh_range = new_state.grid_f->Max() + 1.0; + } + else + { + new_state.SetQuadSolution(); + } } if (curr_state.SetNewMeshAndSolution(std::move(new_state), *vs)) { @@ -859,6 +884,85 @@ void communication_thread::execute() goto comm_terminate; } } + else if (ident == "quadrature" || ident == "pquadrature") + { + bool fix_elem_orient = glvis_command->FixElementOrientations(); + StreamState tmp; + if (ident == "quadrature") + { + tmp.mesh.reset(new Mesh(*is[0], 1, 0, fix_elem_orient)); + tmp.mesh_quad.reset(); + if (!(*is[0])) + { + break; + } + tmp.quad_f.reset(new QuadratureFunction(tmp.mesh.get(), *is[0])); + if (!(*is[0])) + { + break; + } + } + else if (ident == "pquadrature") + { + Array mesh_array; + Array qf_array; + int proc, nproc, np = 0; + bool keep_attr = glvis_command->KeepAttrib(); + do + { + istream &isock = *is[np]; + isock >> nproc >> proc >> ws; +#ifdef GLVIS_DEBUG + cout << "connection[" << np << "]: pquadrature " << nproc << ' ' + << proc << endl; +#endif + isock >> ident >> ws; // "quadrature" + mesh_array.SetSize(nproc); + qf_array.SetSize(nproc); + mesh_array[proc] = new Mesh(isock, 1, 0, fix_elem_orient); + if (!keep_attr) + { + // set element and boundary attributes to proc+1 + for (int i = 0; i < mesh_array[proc]->GetNE(); i++) + { + mesh_array[proc]->GetElement(i)->SetAttribute(proc+1); + } + for (int i = 0; i < mesh_array[proc]->GetNBE(); i++) + { + mesh_array[proc]->GetBdrElement(i)->SetAttribute(proc+1); + } + } + qf_array[proc] = new QuadratureFunction(mesh_array[proc], isock); + np++; + if (np == nproc) + { + break; + } + *is[np] >> ident >> ws; // "pquadrature" + } + while (1); + tmp.mesh.reset(new Mesh(mesh_array, nproc)); + tmp.CollectQuadratures(qf_array, nproc); + + for (int p = 0; p < nproc; p++) + { + delete qf_array[nproc-1-p]; + delete mesh_array[nproc-1-p]; + } + qf_array.DeleteAll(); + mesh_array.DeleteAll(); + } + + // cout << "Stream: new solution" << endl; + + tmp.Extrude1DMeshAndSolution(); + + if (glvis_command->NewMeshAndQuadrature(std::move(tmp.mesh), + std::move(tmp.quad_f))) + { + goto comm_terminate; + } + } else if (ident == "screenshot") { string filename; diff --git a/lib/threads.hpp b/lib/threads.hpp index 804ef7e8..e69d63e1 100644 --- a/lib/threads.hpp +++ b/lib/threads.hpp @@ -108,6 +108,8 @@ class GLVisCommand // called by worker threads int NewMeshAndSolution(std::unique_ptr _new_m, std::unique_ptr _new_g); + int NewMeshAndQuadrature(std::unique_ptr _new_m, + std::unique_ptr _new_q); int Screenshot(const char *filename); int KeyCommands(const char *keys); int WindowSize(int w, int h); From 9bb5f762a5a638a80f5d260ec16c293fff754561 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 1 Jul 2024 13:46:35 -0700 Subject: [PATCH 21/57] Added state of the quadrature representation. --- lib/stream_reader.cpp | 2 ++ lib/stream_reader.hpp | 9 ++++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index a36e7529..ad7aecf8 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -195,6 +195,8 @@ void StreamState::SetQuadSolution(QuadSolution type) gf->MakeOwner(fec); grid_f.reset(gf); } + + quad_sol = type; } // Read the content of an input stream (e.g. from socket/file) diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index ecaad001..93f7ceea 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -49,9 +49,13 @@ struct StreamState enum class QuadSolution { + MIN = -1, + //---------- LOR_GLL, HO_L2, - }; + //---------- + MAX + } quad_sol; /// Helper function for visualizing 1D data void Extrude1DMeshAndSolution(); @@ -65,6 +69,9 @@ struct StreamState /// Set a quadrature function solution producing a proxy grid function void SetQuadSolution(QuadSolution type = QuadSolution::LOR_GLL); + /// Get the current representation of quadrature solution + inline QuadSolution GetQuadSolution() const { return quad_sol; } + FieldType ReadStream(std::istream &is, const std::string &data_type); FieldType ReadStreams(const StreamCollection& input_streams); From f77f7809ea01e148c109fff6c36f4d259d335008 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 1 Jul 2024 13:59:08 -0700 Subject: [PATCH 22/57] Added a key (Q) for switching representation of quad functions. --- glvis.cpp | 17 +++++++++ lib/stream_reader.cpp | 81 ++++++++++++++++++++++--------------------- lib/stream_reader.hpp | 6 ++++ 3 files changed, 65 insertions(+), 39 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index 0747af46..4c444280 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -113,6 +113,9 @@ int ReadParMeshAndGridFunction(int np, const char *mesh_prefix, int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, const char *sol_prefix, StreamState& state); +// switch representation of the quadrature function +void SwitchQuadSolution(); + // Visualize the data in the global variables mesh, sol/grid_f, etc bool GLVisInitVis(StreamState::FieldType field_type, StreamCollection input_streams) @@ -139,6 +142,11 @@ bool GLVisInitVis(StreamState::FieldType field_type, comm_thread = new communication_thread(std::move(input_streams), glvis_command); } + if (stream_state.quad_f) + { + GetAppWindow()->setOnKeyDown('Q', SwitchQuadSolution); + } + double mesh_range = -1.0; if (field_type == StreamState::FieldType::SCALAR || field_type == StreamState::FieldType::MESH) @@ -1940,3 +1948,12 @@ int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, return read_err; } + +void SwitchQuadSolution() +{ + int iqs = ((int)stream_state.GetQuadSolution()+1) + % ((int)StreamState::QuadSolution::MAX); + stream_state.SetQuadSolution((StreamState::QuadSolution)iqs); + stream_state.ResetMeshAndSolution(vs); + SendExposeEvent(); +} diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index ad7aecf8..abb97967 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -572,53 +572,56 @@ bool StreamState::SetNewMeshAndSolution(StreamState new_state, if (new_state.mesh->SpaceDimension() == mesh->SpaceDimension() && new_state.grid_f->VectorDim() == grid_f->VectorDim()) { - std::unique_ptr new_m = std::move(new_state.mesh); - std::unique_ptr new_g = std::move(new_state.grid_f); - std::unique_ptr new_mq = std::move(new_state.mesh_quad); - std::unique_ptr new_q = std::move(new_state.quad_f); - if (new_m->SpaceDimension() == 2) + grid_f = std::move(new_state.grid_f); + mesh = std::move(new_state.mesh); + quad_f = std::move(new_state.quad_f); + mesh_quad = std::move(new_state.mesh_quad); + + ResetMeshAndSolution(vs); + + return true; + } + else + { + return false; + } +} + +void StreamState::ResetMeshAndSolution(VisualizationScene* vs) +{ + if (mesh->SpaceDimension() == 2) + { + if (grid_f->VectorDim() == 1) { - if (new_g->VectorDim() == 1) - { - VisualizationSceneSolution *vss = - dynamic_cast(vs); - new_g->GetNodalValues(sol); - vss->NewMeshAndSolution(new_m.get(), &sol, new_g.get()); - } - else - { - VisualizationSceneVector *vsv = - dynamic_cast(vs); - vsv->NewMeshAndSolution(*new_g); - } + VisualizationSceneSolution *vss = + dynamic_cast(vs); + grid_f->GetNodalValues(sol); + vss->NewMeshAndSolution(mesh.get(), &sol, grid_f.get()); } else { - if (new_g->VectorDim() == 1) - { - VisualizationSceneSolution3d *vss = - dynamic_cast(vs); - new_g->GetNodalValues(sol); - vss->NewMeshAndSolution(new_m.get(), &sol, new_g.get()); - } - else - { - new_g = ProjectVectorFEGridFunction(std::move(new_g)); - - VisualizationSceneVector3d *vss = - dynamic_cast(vs); - vss->NewMeshAndSolution(new_m.get(), new_g.get()); - } + VisualizationSceneVector *vsv = + dynamic_cast(vs); + vsv->NewMeshAndSolution(*grid_f); } - grid_f = std::move(new_g); - mesh = std::move(new_m); - quad_f = std::move(new_q); - mesh_quad = std::move(new_mq); - return true; } else { - return false; + if (grid_f->VectorDim() == 1) + { + VisualizationSceneSolution3d *vss = + dynamic_cast(vs); + grid_f->GetNodalValues(sol); + vss->NewMeshAndSolution(mesh.get(), &sol, grid_f.get()); + } + else + { + grid_f = ProjectVectorFEGridFunction(std::move(grid_f)); + + VisualizationSceneVector3d *vss = + dynamic_cast(vs); + vss->NewMeshAndSolution(mesh.get(), grid_f.get()); + } } } diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index 93f7ceea..b3c241ce 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -87,6 +87,12 @@ struct StreamState /// updated. bool SetNewMeshAndSolution(StreamState new_state, VisualizationScene* vs); + + /// Updates the given VisualizationScene pointer with the new data + /// of this object. + /// @note: Use with caution when the update is compatible + /// @see SetNewMeshAndSolution() + void ResetMeshAndSolution(VisualizationScene* vs); }; // Replace a given VectorFiniteElement-based grid function (e.g. from a Nedelec From e6b0e80deb6ae926152ceb6a09027c7600238b5d Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 1 Jul 2024 14:35:16 -0700 Subject: [PATCH 23/57] Fixed extrusion of 1D quads. --- lib/stream_reader.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index abb97967..ff53f6e0 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -28,6 +28,12 @@ void StreamState::Extrude1DMeshAndSolution() return; } + if (mesh_quad) + { + mesh.swap(mesh_quad); + mesh_quad.reset(); + } + // find xmin and xmax over the vertices of the 1D mesh double xmin = numeric_limits::infinity(); double xmax = -xmin; @@ -52,7 +58,6 @@ void StreamState::Extrude1DMeshAndSolution() Extrude1DQuadFunction(mesh.get(), mesh2d, quad_f.get(), 1); quad_f.reset(quad_f_2d); - SetQuadSolution(); } else if (grid_f) { @@ -72,6 +77,7 @@ void StreamState::Extrude1DMeshAndSolution() } mesh.reset(mesh2d); + if (quad_f) { SetQuadSolution(); } } void StreamState::CollectQuadratures(QuadratureFunction *qf_array[], From 94ab3c186294333c6f37ad5b744324350364aa28 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 1 Jul 2024 15:05:48 -0700 Subject: [PATCH 24/57] Changed extrusion of quads to extrusion of its representation. --- glvis.cpp | 1 + lib/stream_reader.cpp | 17 ++--------------- lib/threads.cpp | 1 + 3 files changed, 4 insertions(+), 15 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index 4c444280..550864a0 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -1954,6 +1954,7 @@ void SwitchQuadSolution() int iqs = ((int)stream_state.GetQuadSolution()+1) % ((int)StreamState::QuadSolution::MAX); stream_state.SetQuadSolution((StreamState::QuadSolution)iqs); + stream_state.Extrude1DMeshAndSolution(); stream_state.ResetMeshAndSolution(vs); SendExposeEvent(); } diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index ff53f6e0..70319a63 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -28,12 +28,6 @@ void StreamState::Extrude1DMeshAndSolution() return; } - if (mesh_quad) - { - mesh.swap(mesh_quad); - mesh_quad.reset(); - } - // find xmin and xmax over the vertices of the 1D mesh double xmin = numeric_limits::infinity(); double xmax = -xmin; @@ -52,14 +46,7 @@ void StreamState::Extrude1DMeshAndSolution() Mesh *mesh2d = Extrude1D(mesh.get(), 1, 0.1*(xmax - xmin)); - if (quad_f) - { - QuadratureFunction *quad_f_2d = - Extrude1DQuadFunction(mesh.get(), mesh2d, quad_f.get(), 1); - - quad_f.reset(quad_f_2d); - } - else if (grid_f) + if (grid_f) { GridFunction *grid_f_2d = Extrude1DGridFunction(mesh.get(), mesh2d, grid_f.get(), 1); @@ -76,8 +63,8 @@ void StreamState::Extrude1DMeshAndSolution() sol = sol2d; } + if (!mesh_quad) { mesh.swap(mesh_quad); } mesh.reset(mesh2d); - if (quad_f) { SetQuadSolution(); } } void StreamState::CollectQuadratures(QuadratureFunction *qf_array[], diff --git a/lib/threads.cpp b/lib/threads.cpp index 4697ce79..54bb9e02 100644 --- a/lib/threads.cpp +++ b/lib/threads.cpp @@ -454,6 +454,7 @@ int GLVisCommand::Execute() else { new_state.SetQuadSolution(); + new_state.Extrude1DMeshAndSolution(); } } if (curr_state.SetNewMeshAndSolution(std::move(new_state), *vs)) From 0a095cb34d0b393b86d24bc2b1fab304125ce9e1 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 1 Jul 2024 16:00:49 -0700 Subject: [PATCH 25/57] Switched LOR refinement to ClosedGL. --- lib/stream_reader.cpp | 4 ++-- lib/stream_reader.hpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 70319a63..56870c81 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -151,13 +151,13 @@ void StreamState::SetQuadSolution(QuadSolution type) mesh.swap(mesh_quad); mesh_quad.reset(); } - if (type == QuadSolution::LOR_GLL) + if (type == QuadSolution::LOR_ClosedGL) { const int ref_factor = order + 1; if (ref_factor > 1) { Mesh *mesh_lor = new Mesh(Mesh::MakeRefined(*mesh, ref_factor, - BasisType::GaussLobatto)); + BasisType::ClosedGL)); FiniteElementCollection *fec = new L2_FECollection(0, mesh->Dimension()); FiniteElementSpace *fes = new FiniteElementSpace(mesh_lor, fec, quad_f->GetVDim(), Ordering::byVDIM); diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index b3c241ce..3c0a4a7f 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -51,7 +51,7 @@ struct StreamState { MIN = -1, //---------- - LOR_GLL, + LOR_ClosedGL, HO_L2, //---------- MAX @@ -67,7 +67,7 @@ struct StreamState void SetMeshSolution(); /// Set a quadrature function solution producing a proxy grid function - void SetQuadSolution(QuadSolution type = QuadSolution::LOR_GLL); + void SetQuadSolution(QuadSolution type = QuadSolution::LOR_ClosedGL); /// Get the current representation of quadrature solution inline QuadSolution GetQuadSolution() const { return quad_sol; } From 427a7621af01a851eba2eae120f3451e6807cd40 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 1 Jul 2024 16:02:17 -0700 Subject: [PATCH 26/57] Added a note to README and CHANGELOG. --- CHANGELOG | 7 +++++++ README.md | 1 + 2 files changed, 8 insertions(+) diff --git a/CHANGELOG b/CHANGELOG index 60cedb5f..d827fbe1 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -21,6 +21,13 @@ Version 4.2.1 (development) - Added a compilation parameter for the default font size. +- Added visualization of quadratures. They can be loaded from a file through the + new command line argument '-q' or in a socket stream with the keyword 'quadrature' + (instead of solution). Two options of visualization are offered (which can be + switched by 'Q' key): piece-wise constant function on a refined mesh (closed + Gauss-Legendre refinement), or discontinuous elements with DOFs collocated with + the quadrature on the original mesh. + Version 4.2 released on May 23, 2022 ==================================== diff --git a/README.md b/README.md index f194ebcf..3c6c806f 100644 --- a/README.md +++ b/README.md @@ -144,6 +144,7 @@ Key commands - Ctrl + p – Print to a PDF file using `gl2ps`. Other vector formats (SVG, EPS) are also possible, but keep in mind that the printing takes a while and the generated files are big. +- Q – Cycle between representations of the quadrature (piece-wise constant refined/L2 element dof collocation) - q – Exit ## Advanced From 6fa2317c0ca2168e03ae774fb95b37e799d31bb1 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 1 Jul 2024 16:29:44 -0700 Subject: [PATCH 27/57] Added a check for tensor-product bases with HO representation. --- lib/stream_reader.cpp | 55 +++++++++++++++++++++++++++++-------------- 1 file changed, 37 insertions(+), 18 deletions(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 56870c81..8c531e5d 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -151,11 +151,17 @@ void StreamState::SetQuadSolution(QuadSolution type) mesh.swap(mesh_quad); mesh_quad.reset(); } - if (type == QuadSolution::LOR_ClosedGL) + + switch (type) { - const int ref_factor = order + 1; - if (ref_factor > 1) + case QuadSolution::LOR_ClosedGL: { + const int ref_factor = order + 1; + if (ref_factor <= 1) + { + SetQuadSolution(QuadSolution::HO_L2);//low-order + return; + } Mesh *mesh_lor = new Mesh(Mesh::MakeRefined(*mesh, ref_factor, BasisType::ClosedGL)); FiniteElementCollection *fec = new L2_FECollection(0, mesh->Dimension()); @@ -170,23 +176,36 @@ void StreamState::SetQuadSolution(QuadSolution type) mesh.swap(mesh_quad); mesh.reset(mesh_lor); } - else + break; + case QuadSolution::HO_L2: { - //low-order - type = QuadSolution::HO_L2; + FiniteElementCollection *fec = new L2_FECollection(order, mesh->Dimension()); + if (order > 0) + { + Array geoms; + mesh->GetGeometries(mesh->Dimension(), geoms); + for (auto geom : geoms) + if (!Geometry::IsTensorProduct(geom)) + { + cout << "High-order representation is available only for tensor-product bases" + << endl; + SetQuadSolution(QuadSolution::LOR_ClosedGL); + return; + } + } + FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, + quad_f->GetVDim(), Ordering::byVDIM); + MFEM_ASSERT(quad_f->Size() == fes->GetVSize(), "Size mismatch"); + cout << "Representing quadrature by grid function of order " + << order << endl; + GridFunction *gf = new GridFunction(fes, *quad_f, 0); + gf->MakeOwner(fec); + grid_f.reset(gf); } - } - if (type == QuadSolution::HO_L2) - { - FiniteElementCollection *fec = new L2_FECollection(order, mesh->Dimension()); - FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, - quad_f->GetVDim(), Ordering::byVDIM); - MFEM_ASSERT(quad_f->Size() == fes->GetVSize(), "Size mismatch"); - cout << "Representing quadrature by grid function of order " - << order << endl; - GridFunction *gf = new GridFunction(fes, *quad_f, 0); - gf->MakeOwner(fec); - grid_f.reset(gf); + break; + default: + cout << "Unknon quadrarture function representation" << endl; + return; } quad_sol = type; From 0b9c3320620f0c2fc4961d563a65af0da129406c Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 2 Jul 2024 14:09:27 -0700 Subject: [PATCH 28/57] Added default quadrature solution type. --- lib/stream_reader.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index 3c0a4a7f..9756aff9 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -49,13 +49,14 @@ struct StreamState enum class QuadSolution { + NONE = -1, MIN = -1, //---------- LOR_ClosedGL, HO_L2, //---------- MAX - } quad_sol; + } quad_sol {QuadSolution::NONE}; /// Helper function for visualizing 1D data void Extrude1DMeshAndSolution(); From d52af7f1d61e780c852b82dffce85ef207d7ae75 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 2 Jul 2024 15:59:55 -0700 Subject: [PATCH 29/57] Protected writes to the pointer members of StreamState to keep consistency. --- glvis.cpp | 53 +++++++-------- lib/stream_reader.cpp | 147 +++++++++++++++++++++++++++++++----------- lib/stream_reader.hpp | 46 ++++++++++--- lib/threads.cpp | 49 ++++---------- lib/threads.hpp | 5 +- 5 files changed, 182 insertions(+), 118 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index 550864a0..056abbb1 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -238,8 +238,7 @@ bool GLVisInitVis(StreamState::FieldType field_type, { if (stream_state.grid_f) { - stream_state.grid_f - = ProjectVectorFEGridFunction(std::move(stream_state.grid_f)); + stream_state.ProjectVectorFEGridFunction(); vs = new VisualizationSceneVector3d(*stream_state.grid_f); } else @@ -304,14 +303,14 @@ int ScriptReadSolution(istream &scr, StreamState& state) cout << "Can not open mesh file: " << mword << endl; return 1; } - state.mesh.reset(new Mesh(imesh, 1, 0, state.fix_elem_orient)); + state.SetMesh(new Mesh(imesh, 1, 0, state.fix_elem_orient)); // read the solution (GridFunction) scr >> ws >> sword; if (sword == mword) // mesh and solution in the same file { cout << "solution: " << mword << endl; - state.grid_f.reset(new GridFunction(state.mesh.get(), imesh)); + state.SetGridFunction(new GridFunction(state.mesh.get(), imesh)); } else { @@ -322,7 +321,7 @@ int ScriptReadSolution(istream &scr, StreamState& state) cout << "Can not open solution file: " << sword << endl; return 2; } - state.grid_f.reset(new GridFunction(state.mesh.get(), isol)); + state.SetGridFunction(new GridFunction(state.mesh.get(), isol)); } state.Extrude1DMeshAndSolution(); @@ -344,15 +343,14 @@ int ScriptReadQuadrature(istream &scr, StreamState& state) cout << "Can not open mesh file: " << mword << endl; return 1; } - state.mesh.reset(new Mesh(imesh, 1, 0, state.fix_elem_orient)); - state.mesh_quad.reset(); + state.SetMesh(new Mesh(imesh, 1, 0, state.fix_elem_orient)); // read the quadrature (QuadratureFunction) scr >> ws >> sword; if (sword == mword) // mesh and quadrature in the same file { cout << "quadrature: " << mword << endl; - state.quad_f.reset(new QuadratureFunction(state.mesh.get(), imesh)); + state.SetQuadFunction(new QuadratureFunction(state.mesh.get(), imesh)); } else { @@ -363,7 +361,7 @@ int ScriptReadQuadrature(istream &scr, StreamState& state) cout << "Can not open quadrature file: " << sword << endl; return 2; } - state.quad_f.reset(new QuadratureFunction(state.mesh.get(), isol)); + state.SetQuadFunction(new QuadratureFunction(state.mesh.get(), isol)); } state.Extrude1DMeshAndSolution(); @@ -454,7 +452,7 @@ int ScriptReadDisplMesh(istream &scr, StreamState& state) return 1; } cout << word << endl; - meshstate.mesh.reset(new Mesh(imesh, 1, 0, state.fix_elem_orient)); + meshstate.SetMesh(new Mesh(imesh, 1, 0, state.fix_elem_orient)); } meshstate.Extrude1DMeshAndSolution(); Mesh* const m = meshstate.mesh.get(); @@ -462,8 +460,8 @@ int ScriptReadDisplMesh(istream &scr, StreamState& state) { init_nodes = new Vector; meshstate.mesh->GetNodes(*init_nodes); - state.mesh = NULL; - state.grid_f = NULL; + state.SetMesh(NULL); + state.SetGridFunction(NULL); } else { @@ -476,7 +474,7 @@ int ScriptReadDisplMesh(istream &scr, StreamState& state) vfes = new FiniteElementSpace(m, vfec, m->SpaceDimension()); } - meshstate.grid_f.reset(new GridFunction(vfes)); + meshstate.SetGridFunction(new GridFunction(vfes)); GridFunction * const g = meshstate.grid_f.get(); if (vfec) { @@ -493,8 +491,7 @@ int ScriptReadDisplMesh(istream &scr, StreamState& state) *g = 0.0; } - state.mesh = std::move(meshstate.mesh); - state.grid_f = std::move(meshstate.grid_f); + state = std::move(meshstate); } return 0; @@ -1598,8 +1595,7 @@ void ReadSerial(StreamState& state) exit(1); } - state.mesh.reset(new Mesh(meshin, 1, 0, state.fix_elem_orient)); - state.mesh_quad.reset(); + state.SetMesh(new Mesh(meshin, 1, 0, state.fix_elem_orient)); if (state.is_gf || state.is_qf || (input & INPUT_SCALAR_SOL) || (input & INPUT_VECTOR_SOL)) @@ -1624,12 +1620,12 @@ void ReadSerial(StreamState& state) if (state.is_gf) { - state.grid_f.reset(new GridFunction(state.mesh.get(), *solin)); + state.SetGridFunction(new GridFunction(state.mesh.get(), *solin)); SetGridFunction(state); } else if (state.is_qf) { - state.quad_f.reset(new QuadratureFunction(state.mesh.get(), *solin)); + state.SetQuadFunction(new QuadratureFunction(state.mesh.get(), *solin)); SetQuadFunction(state); } else if (input & INPUT_SCALAR_SOL) @@ -1681,7 +1677,7 @@ void SetGridFunction(StreamState& state) { (*new_gf)(i) = (*state.grid_f)(ofes->DofToVDof(i, gf_component)); } - state.grid_f.reset(new_gf); + state.SetGridFunction(new_gf); } if (state.grid_f->VectorDim() == 1) { @@ -1712,7 +1708,7 @@ void SetQuadFunction(StreamState& state) } state.quad_f->SetOwnsSpace(false); new_qf->SetOwnsSpace(true); - state.quad_f.reset(new_qf); + state.SetQuadFunction(new_qf); } if (vdim == 1) { @@ -1770,14 +1766,14 @@ int ReadParMeshAndGridFunction(int np, const char *mesh_prefix, const char *sol_prefix, StreamState& state) { - state.mesh = NULL; + state.SetMesh(NULL); // are the solutions bundled together with the mesh files? bool same_file = false; if (sol_prefix) { same_file = !strcmp(sol_prefix, mesh_prefix); - state.grid_f = NULL; + state.SetGridFunction(NULL); } Array mesh_array(np); @@ -1841,10 +1837,10 @@ int ReadParMeshAndGridFunction(int np, const char *mesh_prefix, if (!read_err) { // create the combined mesh and gf - state.mesh.reset(new Mesh(mesh_array, np)); + state.SetMesh(new Mesh(mesh_array, np)); if (sol_prefix) { - state.grid_f.reset(new GridFunction(state.mesh.get(), gf_array, np)); + state.SetGridFunction(new GridFunction(state.mesh.get(), gf_array, np)); } } @@ -1861,14 +1857,14 @@ int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, const char *sol_prefix, StreamState& state) { - state.mesh = NULL; + state.SetMesh(NULL); // are the solutions bundled together with the mesh files? bool same_file = false; if (sol_prefix) { same_file = !strcmp(sol_prefix, mesh_prefix); - state.grid_f = NULL; + state.SetGridFunction(NULL); } Array mesh_array(np); @@ -1932,8 +1928,7 @@ int ReadParMeshAndQuadFunction(int np, const char *mesh_prefix, if (!read_err) { // create the combined mesh and gf - state.mesh.reset(new Mesh(mesh_array, np)); - state.mesh_quad.reset(); + state.SetMesh(new Mesh(mesh_array, np)); if (sol_prefix) { state.CollectQuadratures(qf_array, np); diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 8c531e5d..12d3e969 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -21,6 +21,77 @@ using namespace mfem; QuadratureFunction* Extrude1DQuadFunction(Mesh *mesh, Mesh *mesh2d, QuadratureFunction *qf, int ny); +StreamState &StreamState::operator=(StreamState &&ss) +{ + internal.mesh = move(ss.internal.mesh); + internal.mesh_quad = move(ss.internal.mesh_quad); + internal.grid_f = move(ss.internal.grid_f); + internal.quad_f = move(ss.internal.quad_f); + + sol = move(ss.sol); + solu = move(ss.solu); + solv = move(ss.solv); + solw = move(ss.solw); + normals = move(ss.normals); + keys = move(ss.keys); + + is_gf = ss.is_gf; + is_qf = ss.is_qf; + fix_elem_orient = ss.fix_elem_orient; + save_coloring = ss.save_coloring; + keep_attr = ss.keep_attr; + + return *this; +} + +void StreamState::SetMesh(mfem::Mesh *mesh_) +{ + internal.mesh.reset(mesh_); + SetMesh(move(internal.mesh)); +} + +void StreamState::SetMesh(std::unique_ptr &&pmesh) +{ + internal.mesh = move(pmesh); + internal.mesh_quad.reset(); + if (grid_f && grid_f->FESpace()->GetMesh() != mesh.get()) { SetGridFunction(NULL); } + if (quad_f && quad_f->GetSpace()->GetMesh() != mesh.get()) { SetQuadFunction(NULL); } +} + +void StreamState::SetGridFunction(mfem::GridFunction *gf) +{ + internal.grid_f.reset(gf); + SetGridFunction(move(internal.grid_f)); +} + +void StreamState::SetGridFunction(std::unique_ptr &&pgf) +{ + internal.grid_f = move(pgf); + internal.quad_f.reset(); + quad_sol = QuadSolution::NONE; +} + +void StreamState::SetQuadFunction(mfem::QuadratureFunction *qf) +{ + if (quad_f.get() != qf) + { + internal.grid_f.reset(); + quad_sol = QuadSolution::NONE; + } + internal.quad_f.reset(qf); +} + +void StreamState::SetQuadFunction(std::unique_ptr + &&pqf) +{ + if (quad_f.get() != pqf.get()) + { + internal.grid_f.reset(); + quad_sol = QuadSolution::NONE; + } + internal.quad_f = move(pqf); +} + void StreamState::Extrude1DMeshAndSolution() { if (mesh->Dimension() != 1 || mesh->SpaceDimension() != 1) @@ -51,7 +122,7 @@ void StreamState::Extrude1DMeshAndSolution() GridFunction *grid_f_2d = Extrude1DGridFunction(mesh.get(), mesh2d, grid_f.get(), 1); - grid_f.reset(grid_f_2d); + internal.grid_f.reset(grid_f_2d); } else if (sol.Size() == mesh->GetNV()) { @@ -63,8 +134,8 @@ void StreamState::Extrude1DMeshAndSolution() sol = sol2d; } - if (!mesh_quad) { mesh.swap(mesh_quad); } - mesh.reset(mesh2d); + if (!mesh_quad) { internal.mesh.swap(internal.mesh_quad); } + internal.mesh.reset(mesh2d); } void StreamState::CollectQuadratures(QuadratureFunction *qf_array[], @@ -75,7 +146,7 @@ void StreamState::CollectQuadratures(QuadratureFunction *qf_array[], //assume the same quadrature rule QuadratureSpace *qspace = new QuadratureSpace(*mesh, qf_array[0]->GetIntRule(0)); - quad_f.reset(new QuadratureFunction(qspace, vdim)); + SetQuadFunction(new QuadratureFunction(qspace, vdim)); quad_f->SetOwnsSpace(true); real_t *g_data = quad_f->GetData(); for (int p = 0; p < npieces; p++) @@ -107,7 +178,7 @@ void StreamState::SetMeshSolution() cfec = new Const3DFECollection; } FiniteElementSpace *cfes = new FiniteElementSpace(mesh.get(), cfec); - grid_f.reset(new GridFunction(cfes)); + SetGridFunction(new GridFunction(cfes)); grid_f->MakeOwner(cfec); { Array coloring; @@ -148,8 +219,8 @@ void StreamState::SetQuadSolution(QuadSolution type) //use the original mesh when available if (mesh_quad.get()) { - mesh.swap(mesh_quad); - mesh_quad.reset(); + internal.mesh.swap(internal.mesh_quad); + internal.mesh_quad.reset(); } switch (type) @@ -172,9 +243,9 @@ void StreamState::SetQuadSolution(QuadSolution type) << ref_factor << " times" << endl; GridFunction *gf = new GridFunction(fes, *quad_f, 0); gf->MakeOwner(fec); - grid_f.reset(gf); - mesh.swap(mesh_quad); - mesh.reset(mesh_lor); + internal.grid_f.reset(gf); + internal.mesh.swap(internal.mesh_quad); + internal.mesh.reset(mesh_lor); } break; case QuadSolution::HO_L2: @@ -200,7 +271,7 @@ void StreamState::SetQuadSolution(QuadSolution type) << order << endl; GridFunction *gf = new GridFunction(fes, *quad_f, 0); gf->MakeOwner(fec); - grid_f.reset(gf); + internal.grid_f.reset(gf); } break; default: @@ -219,13 +290,13 @@ StreamState::FieldType StreamState::ReadStream(istream &is, keys.clear(); if (data_type == "fem2d_data") { - mesh.reset(new Mesh(is, 0, 0, fix_elem_orient)); + SetMesh(new Mesh(is, 0, 0, fix_elem_orient)); sol.Load(is, mesh->GetNV()); } else if (data_type == "vfem2d_data" || data_type == "vfem2d_data_keys") { field_type = FieldType::VECTOR; - mesh.reset(new Mesh(is, 0, 0, fix_elem_orient)); + SetMesh(new Mesh(is, 0, 0, fix_elem_orient)); solu.Load(is, mesh->GetNV()); solv.Load(is, mesh->GetNV()); if (data_type == "vfem2d_data_keys") @@ -235,13 +306,13 @@ StreamState::FieldType StreamState::ReadStream(istream &is, } else if (data_type == "fem3d_data") { - mesh.reset(new Mesh(is, 0, 0, fix_elem_orient)); + SetMesh(new Mesh(is, 0, 0, fix_elem_orient)); sol.Load(is, mesh->GetNV()); } else if (data_type == "vfem3d_data" || data_type == "vfem3d_data_keys") { field_type = FieldType::VECTOR; - mesh.reset(new Mesh(is, 0, 0, fix_elem_orient)); + SetMesh(new Mesh(is, 0, 0, fix_elem_orient)); solu.Load(is, mesh->GetNV()); solv.Load(is, mesh->GetNV()); solw.Load(is, mesh->GetNV()); @@ -252,8 +323,8 @@ StreamState::FieldType StreamState::ReadStream(istream &is, } else if (data_type == "fem2d_gf_data" || data_type == "fem2d_gf_data_keys") { - mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); - grid_f.reset(new GridFunction(mesh.get(), is)); + SetMesh(new Mesh(is, 1, 0, fix_elem_orient)); + SetGridFunction(new GridFunction(mesh.get(), is)); if (data_type == "fem2d_gf_data_keys") { is >> keys; @@ -262,8 +333,8 @@ StreamState::FieldType StreamState::ReadStream(istream &is, else if (data_type == "vfem2d_gf_data" || data_type == "vfem2d_gf_data_keys") { field_type = FieldType::VECTOR; - mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); - grid_f.reset(new GridFunction(mesh.get(), is)); + SetMesh(new Mesh(is, 1, 0, fix_elem_orient)); + SetGridFunction(new GridFunction(mesh.get(), is)); if (data_type == "vfem2d_gf_data_keys") { is >> keys; @@ -271,8 +342,8 @@ StreamState::FieldType StreamState::ReadStream(istream &is, } else if (data_type == "fem3d_gf_data" || data_type == "fem3d_gf_data_keys") { - mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); - grid_f.reset(new GridFunction(mesh.get(), is)); + SetMesh(new Mesh(is, 1, 0, fix_elem_orient)); + SetGridFunction(new GridFunction(mesh.get(), is)); if (data_type == "fem3d_gf_data_keys") { is >> keys; @@ -281,8 +352,8 @@ StreamState::FieldType StreamState::ReadStream(istream &is, else if (data_type == "vfem3d_gf_data" || data_type == "vfem3d_gf_data_keys") { field_type = FieldType::VECTOR; - mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); - grid_f.reset(new GridFunction(mesh.get(), is)); + SetMesh(new Mesh(is, 1, 0, fix_elem_orient)); + SetGridFunction(new GridFunction(mesh.get(), is)); if (data_type == "vfem3d_gf_data_keys") { is >> keys; @@ -290,21 +361,20 @@ StreamState::FieldType StreamState::ReadStream(istream &is, } else if (data_type == "solution") { - mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); - grid_f.reset(new GridFunction(mesh.get(), is)); + SetMesh(new Mesh(is, 1, 0, fix_elem_orient)); + SetGridFunction(new GridFunction(mesh.get(), is)); field_type = (grid_f->VectorDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR; } else if (data_type == "quadrature") { - mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); - mesh_quad.reset(); - quad_f.reset(new QuadratureFunction(mesh.get(), is)); + SetMesh(new Mesh(is, 1, 0, fix_elem_orient)); + SetQuadFunction(new QuadratureFunction(mesh.get(), is)); SetQuadSolution(); field_type = (quad_f->GetVDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR; } else if (data_type == "mesh") { - mesh.reset(new Mesh(is, 1, 0, fix_elem_orient)); + SetMesh(new Mesh(is, 1, 0, fix_elem_orient)); SetMeshSolution(); field_type = FieldType::MESH; } @@ -363,7 +433,7 @@ StreamState::FieldType StreamState::ReadStream(istream &is, tot_num_elem += num_elem; } - mesh.reset(new Mesh(2, tot_num_vert, tot_num_elem, 0)); + SetMesh(new Mesh(2, tot_num_vert, tot_num_elem, 0)); sol.SetSize(tot_num_vert); normals.SetSize(3*tot_num_vert); @@ -498,11 +568,10 @@ StreamState::FieldType StreamState::ReadStreams(const StreamCollection& mfem_error("Input streams contain a mixture of data types!"); } - mesh.reset(new Mesh(mesh_array, nproc)); - mesh_quad.reset(); + SetMesh(new Mesh(mesh_array, nproc)); if (gf_count > 0) { - grid_f.reset(new GridFunction(mesh.get(), gf_array, nproc)); + SetGridFunction(new GridFunction(mesh.get(), gf_array, nproc)); field_type = (grid_f->VectorDim() == 1) ? FieldType::SCALAR : FieldType::VECTOR; } else if (qf_count > 0) @@ -556,7 +625,7 @@ void StreamState::WriteStream(std::ostream &os) // or Raviart-Thomas space) with a discontinuous piece-wise polynomial Cartesian // product vector grid function of the same order. std::unique_ptr -ProjectVectorFEGridFunction(std::unique_ptr gf) +StreamState::ProjectVectorFEGridFunction(std::unique_ptr gf) { if ((gf->VectorDim() == 3) && (gf->FESpace()->GetVDim() == 1)) { @@ -584,10 +653,10 @@ bool StreamState::SetNewMeshAndSolution(StreamState new_state, if (new_state.mesh->SpaceDimension() == mesh->SpaceDimension() && new_state.grid_f->VectorDim() == grid_f->VectorDim()) { - grid_f = std::move(new_state.grid_f); - mesh = std::move(new_state.mesh); - quad_f = std::move(new_state.quad_f); - mesh_quad = std::move(new_state.mesh_quad); + internal.grid_f = std::move(new_state.internal.grid_f); + internal.mesh = std::move(new_state.internal.mesh); + internal.quad_f = std::move(new_state.internal.quad_f); + internal.mesh_quad = std::move(new_state.internal.mesh_quad); ResetMeshAndSolution(vs); @@ -628,7 +697,7 @@ void StreamState::ResetMeshAndSolution(VisualizationScene* vs) } else { - grid_f = ProjectVectorFEGridFunction(std::move(grid_f)); + SetGridFunction(ProjectVectorFEGridFunction(std::move(internal.grid_f))); VisualizationSceneVector3d *vss = dynamic_cast(vs); diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index 9756aff9..ab02f73c 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -23,12 +23,22 @@ using StreamCollection = std::vector>; struct StreamState { +private: + struct + { + std::unique_ptr mesh; + std::unique_ptr mesh_quad; + std::unique_ptr grid_f; + std::unique_ptr quad_f; + } internal; + +public: mfem::Vector sol, solu, solv, solw, normals; std::string keys; - std::unique_ptr mesh; - std::unique_ptr mesh_quad; - std::unique_ptr grid_f; - std::unique_ptr quad_f; + 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 &quad_f{internal.quad_f}; int is_gf{0}; int is_qf{0}; bool fix_elem_orient{false}; @@ -58,6 +68,19 @@ struct StreamState MAX } quad_sol {QuadSolution::NONE}; + StreamState() = default; + StreamState(StreamState &&ss) { *this = std::move(ss); } + StreamState& operator=(StreamState &&ss); + + void SetMesh(mfem::Mesh *mesh); + void SetMesh(std::unique_ptr &&pmesh); + + void SetGridFunction(mfem::GridFunction *gf); + void SetGridFunction(std::unique_ptr &&pgf); + + void SetQuadFunction(mfem::QuadratureFunction *qf); + void SetQuadFunction(std::unique_ptr &&pqf); + /// Helper function for visualizing 1D data void Extrude1DMeshAndSolution(); @@ -79,6 +102,15 @@ struct StreamState void WriteStream(std::ostream &os); + // 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); + + void ProjectVectorFEGridFunction() + { internal.grid_f = ProjectVectorFEGridFunction(std::move(internal.grid_f)); } + /// Sets a new mesh and solution from another StreamState object, and /// updates the given VisualizationScene pointer with the new data. /// @@ -96,10 +128,4 @@ struct StreamState void ResetMeshAndSolution(VisualizationScene* vs); }; -// 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. -std::unique_ptr -ProjectVectorFEGridFunction(std::unique_ptr gf); - #endif // GLVIS_STREAM_READER_HPP diff --git a/lib/threads.cpp b/lib/threads.cpp index 54bb9e02..0e390c40 100644 --- a/lib/threads.cpp +++ b/lib/threads.cpp @@ -79,34 +79,14 @@ void GLVisCommand::unlock() } } -int GLVisCommand::NewMeshAndSolution(std::unique_ptr _new_m, - std::unique_ptr _new_g) +int GLVisCommand::NewMeshAndSolution(StreamState &&ss) { if (lock() < 0) { return -1; } command = NEW_MESH_AND_SOLUTION; - new_state.mesh = std::move(_new_m); - new_state.grid_f = std::move(_new_g); - if (signal() < 0) - { - return -2; - } - return 0; -} - -int GLVisCommand::NewMeshAndQuadrature(std::unique_ptr _new_m, - std::unique_ptr _new_q) -{ - if (lock() < 0) - { - return -1; - } - command = NEW_MESH_AND_SOLUTION; - new_state.mesh = std::move(_new_m); - new_state.mesh_quad.reset(); - new_state.quad_f = std::move(_new_q); + new_state = move(ss); if (signal() < 0) { return -2; @@ -804,21 +784,21 @@ void communication_thread::execute() StreamState tmp; if (ident == "mesh") { - tmp.mesh.reset(new Mesh(*is[0], 1, 0, fix_elem_orient)); + tmp.SetMesh(new Mesh(*is[0], 1, 0, fix_elem_orient)); if (!(*is[0])) { break; } - tmp.grid_f = NULL; + tmp.SetGridFunction(NULL); } else if (ident == "solution") { - tmp.mesh.reset(new Mesh(*is[0], 1, 0, fix_elem_orient)); + tmp.SetMesh(new Mesh(*is[0], 1, 0, fix_elem_orient)); if (!(*is[0])) { break; } - tmp.grid_f.reset(new GridFunction(tmp.mesh.get(), *is[0])); + tmp.SetGridFunction(new GridFunction(tmp.mesh.get(), *is[0])); if (!(*is[0])) { break; @@ -863,8 +843,8 @@ void communication_thread::execute() *is[np] >> ident >> ws; // "parallel" } while (1); - tmp.mesh.reset(new Mesh(mesh_array, nproc)); - tmp.grid_f.reset(new GridFunction(tmp.mesh.get(), gf_array, nproc)); + tmp.SetMesh(new Mesh(mesh_array, nproc)); + tmp.SetGridFunction(new GridFunction(tmp.mesh.get(), gf_array, nproc)); for (int p = 0; p < nproc; p++) { @@ -879,8 +859,7 @@ void communication_thread::execute() tmp.Extrude1DMeshAndSolution(); - if (glvis_command->NewMeshAndSolution(std::move(tmp.mesh), - std::move(tmp.grid_f))) + if (glvis_command->NewMeshAndSolution(std::move(tmp))) { goto comm_terminate; } @@ -891,13 +870,12 @@ void communication_thread::execute() StreamState tmp; if (ident == "quadrature") { - tmp.mesh.reset(new Mesh(*is[0], 1, 0, fix_elem_orient)); - tmp.mesh_quad.reset(); + tmp.SetMesh(new Mesh(*is[0], 1, 0, fix_elem_orient)); if (!(*is[0])) { break; } - tmp.quad_f.reset(new QuadratureFunction(tmp.mesh.get(), *is[0])); + tmp.SetQuadFunction(new QuadratureFunction(tmp.mesh.get(), *is[0])); if (!(*is[0])) { break; @@ -942,7 +920,7 @@ void communication_thread::execute() *is[np] >> ident >> ws; // "pquadrature" } while (1); - tmp.mesh.reset(new Mesh(mesh_array, nproc)); + tmp.SetMesh(new Mesh(mesh_array, nproc)); tmp.CollectQuadratures(qf_array, nproc); for (int p = 0; p < nproc; p++) @@ -958,8 +936,7 @@ void communication_thread::execute() tmp.Extrude1DMeshAndSolution(); - if (glvis_command->NewMeshAndQuadrature(std::move(tmp.mesh), - std::move(tmp.quad_f))) + if (glvis_command->NewMeshAndSolution(std::move(tmp))) { goto comm_terminate; } diff --git a/lib/threads.hpp b/lib/threads.hpp index e69d63e1..8b5bce2e 100644 --- a/lib/threads.hpp +++ b/lib/threads.hpp @@ -106,10 +106,7 @@ class GLVisCommand bool FixElementOrientations() { return curr_state.fix_elem_orient; } // called by worker threads - int NewMeshAndSolution(std::unique_ptr _new_m, - std::unique_ptr _new_g); - int NewMeshAndQuadrature(std::unique_ptr _new_m, - std::unique_ptr _new_q); + int NewMeshAndSolution(StreamState &&ss); int Screenshot(const char *filename); int KeyCommands(const char *keys); int WindowSize(int w, int h); From c63a983feb1036af8af4f8cae16231d6612d5efa Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 2 Jul 2024 16:08:14 -0700 Subject: [PATCH 30/57] Fixed move -> std::move. --- lib/stream_reader.cpp | 32 ++++++++++++++++---------------- lib/threads.cpp | 2 +- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 12d3e969..47278094 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -23,17 +23,17 @@ QuadratureFunction* Extrude1DQuadFunction(Mesh *mesh, Mesh *mesh2d, StreamState &StreamState::operator=(StreamState &&ss) { - internal.mesh = move(ss.internal.mesh); - internal.mesh_quad = move(ss.internal.mesh_quad); - internal.grid_f = move(ss.internal.grid_f); - internal.quad_f = move(ss.internal.quad_f); - - sol = move(ss.sol); - solu = move(ss.solu); - solv = move(ss.solv); - solw = move(ss.solw); - normals = move(ss.normals); - keys = move(ss.keys); + internal.mesh = std::move(ss.internal.mesh); + internal.mesh_quad = std::move(ss.internal.mesh_quad); + internal.grid_f = std::move(ss.internal.grid_f); + internal.quad_f = std::move(ss.internal.quad_f); + + sol = std::move(ss.sol); + solu = std::move(ss.solu); + solv = std::move(ss.solv); + solw = std::move(ss.solw); + normals = std::move(ss.normals); + keys = std::move(ss.keys); is_gf = ss.is_gf; is_qf = ss.is_qf; @@ -47,12 +47,12 @@ StreamState &StreamState::operator=(StreamState &&ss) void StreamState::SetMesh(mfem::Mesh *mesh_) { internal.mesh.reset(mesh_); - SetMesh(move(internal.mesh)); + SetMesh(std::move(internal.mesh)); } void StreamState::SetMesh(std::unique_ptr &&pmesh) { - internal.mesh = move(pmesh); + internal.mesh = std::move(pmesh); internal.mesh_quad.reset(); if (grid_f && grid_f->FESpace()->GetMesh() != mesh.get()) { SetGridFunction(NULL); } if (quad_f && quad_f->GetSpace()->GetMesh() != mesh.get()) { SetQuadFunction(NULL); } @@ -61,12 +61,12 @@ void StreamState::SetMesh(std::unique_ptr &&pmesh) void StreamState::SetGridFunction(mfem::GridFunction *gf) { internal.grid_f.reset(gf); - SetGridFunction(move(internal.grid_f)); + SetGridFunction(std::move(internal.grid_f)); } void StreamState::SetGridFunction(std::unique_ptr &&pgf) { - internal.grid_f = move(pgf); + internal.grid_f = std::move(pgf); internal.quad_f.reset(); quad_sol = QuadSolution::NONE; } @@ -89,7 +89,7 @@ void StreamState::SetQuadFunction(std::unique_ptr internal.grid_f.reset(); quad_sol = QuadSolution::NONE; } - internal.quad_f = move(pqf); + internal.quad_f = std::move(pqf); } void StreamState::Extrude1DMeshAndSolution() diff --git a/lib/threads.cpp b/lib/threads.cpp index 0e390c40..f9e9d512 100644 --- a/lib/threads.cpp +++ b/lib/threads.cpp @@ -86,7 +86,7 @@ int GLVisCommand::NewMeshAndSolution(StreamState &&ss) return -1; } command = NEW_MESH_AND_SOLUTION; - new_state = move(ss); + new_state = std::move(ss); if (signal() < 0) { return -2; From ef8f2946f2f7837331f71b8992c6ac2090718f42 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 8 Jul 2024 14:10:24 -0700 Subject: [PATCH 31/57] Added continuation of the quad representation. --- lib/threads.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/lib/threads.cpp b/lib/threads.cpp index f9e9d512..23d309fe 100644 --- a/lib/threads.cpp +++ b/lib/threads.cpp @@ -433,7 +433,15 @@ int GLVisCommand::Execute() } else { - new_state.SetQuadSolution(); + auto qs = curr_state.GetQuadSolution(); + if (qs != StreamState::QuadSolution::NONE) + { + new_state.SetQuadSolution(qs); + } + else + { + new_state.SetQuadSolution(); + } new_state.Extrude1DMeshAndSolution(); } } From 43ecc5978dd0f152e6a2281a88f2f3349f45f864 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 8 Jul 2024 14:34:52 -0700 Subject: [PATCH 32/57] Fixed and simplified StreamState assignment operator. --- lib/stream_reader.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 47278094..1ecc9a37 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -23,10 +23,7 @@ QuadratureFunction* Extrude1DQuadFunction(Mesh *mesh, Mesh *mesh2d, StreamState &StreamState::operator=(StreamState &&ss) { - internal.mesh = std::move(ss.internal.mesh); - internal.mesh_quad = std::move(ss.internal.mesh_quad); - internal.grid_f = std::move(ss.internal.grid_f); - internal.quad_f = std::move(ss.internal.quad_f); + internal = std::move(ss.internal); sol = std::move(ss.sol); solu = std::move(ss.solu); @@ -41,6 +38,8 @@ StreamState &StreamState::operator=(StreamState &&ss) save_coloring = ss.save_coloring; keep_attr = ss.keep_attr; + quad_sol = ss.quad_sol; + return *this; } From f5451f3bbdade6cdaae018738ef9bc0f506d4581 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 8 Jul 2024 14:56:14 -0700 Subject: [PATCH 33/57] Mentioned tensor basis in changelog. --- CHANGELOG | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG b/CHANGELOG index d827fbe1..beb526db 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -26,7 +26,7 @@ Version 4.2.1 (development) (instead of solution). Two options of visualization are offered (which can be switched by 'Q' key): piece-wise constant function on a refined mesh (closed Gauss-Legendre refinement), or discontinuous elements with DOFs collocated with - the quadrature on the original mesh. + the quadrature on the original mesh (available for tensor basis only). Version 4.2 released on May 23, 2022 From 9ed2a78224cb4d9806e0b57440ed79d169b58536 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 8 Jul 2024 16:31:55 -0700 Subject: [PATCH 34/57] Fixed loading of scripts. --- glvis.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/glvis.cpp b/glvis.cpp index 056abbb1..f9d41508 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -364,6 +364,7 @@ int ScriptReadQuadrature(istream &scr, StreamState& state) state.SetQuadFunction(new QuadratureFunction(state.mesh.get(), isol)); } + state.SetQuadSolution(); state.Extrude1DMeshAndSolution(); return 0; @@ -432,6 +433,7 @@ int ScriptReadParQuadrature(istream &scr, StreamState& state) quad_prefix.c_str(), state); if (!err_read) { + state.SetQuadSolution(); state.Extrude1DMeshAndSolution(); } return err_read; From 2bbee2b741af880acd271a61dfc852cafe4110cf Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 8 Jul 2024 17:24:34 -0700 Subject: [PATCH 35/57] Fixed parallel assert in CollectQuadratures(). --- 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 1ecc9a37..03134f14 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -152,7 +152,7 @@ void StreamState::CollectQuadratures(QuadratureFunction *qf_array[], { const real_t *l_data = qf_array[p]->GetData(); const int l_size = qf_array[p]->Size(); - MFEM_ASSERT(g_data + l_size >= quad_f->GetData() + quad_f->Size(), + MFEM_ASSERT(g_data + l_size <= quad_f->GetData() + quad_f->Size(), "Local parts do not fit to the global quadrature function!"); memcpy(g_data, l_data, l_size * sizeof(real_t)); g_data += l_size; From b1dc68bea77f2418f8d20a01e8fcca679c7eadf6 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 8 Jul 2024 17:50:55 -0700 Subject: [PATCH 36/57] Fixed 3D vector fields projection for quads. --- 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 03134f14..9354c698 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -696,7 +696,7 @@ void StreamState::ResetMeshAndSolution(VisualizationScene* vs) } else { - SetGridFunction(ProjectVectorFEGridFunction(std::move(internal.grid_f))); + ProjectVectorFEGridFunction(); VisualizationSceneVector3d *vss = dynamic_cast(vs); From 8a8c3db7929866ca368be856b8d34a57dff16371 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 9 Jul 2024 09:14:03 -0700 Subject: [PATCH 37/57] Extended the check for tensor-product elements. --- CHANGELOG | 3 ++- lib/stream_reader.cpp | 27 +++++++++++++-------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index beb526db..968b0dbc 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -26,7 +26,8 @@ Version 4.2.1 (development) (instead of solution). Two options of visualization are offered (which can be switched by 'Q' key): piece-wise constant function on a refined mesh (closed Gauss-Legendre refinement), or discontinuous elements with DOFs collocated with - the quadrature on the original mesh (available for tensor basis only). + the quadrature on the original mesh. High-order quadratures are supported only + for tensor finite elements. Version 4.2 released on May 23, 2022 diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 9354c698..8118d2ab 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -222,6 +222,18 @@ void StreamState::SetQuadSolution(QuadSolution type) internal.mesh_quad.reset(); } + //check for tensor-product basis + if (order > 0) + { + Array geoms; + mesh->GetGeometries(mesh->Dimension(), geoms); + for (auto geom : geoms) + if (!Geometry::IsTensorProduct(geom)) + { + MFEM_ABORT("High-order quadratures are available only for tensor-product finite elements"); + } + } + switch (type) { case QuadSolution::LOR_ClosedGL: @@ -250,19 +262,6 @@ void StreamState::SetQuadSolution(QuadSolution type) case QuadSolution::HO_L2: { FiniteElementCollection *fec = new L2_FECollection(order, mesh->Dimension()); - if (order > 0) - { - Array geoms; - mesh->GetGeometries(mesh->Dimension(), geoms); - for (auto geom : geoms) - if (!Geometry::IsTensorProduct(geom)) - { - cout << "High-order representation is available only for tensor-product bases" - << endl; - SetQuadSolution(QuadSolution::LOR_ClosedGL); - return; - } - } FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, quad_f->GetVDim(), Ordering::byVDIM); MFEM_ASSERT(quad_f->Size() == fes->GetVSize(), "Size mismatch"); @@ -274,7 +273,7 @@ void StreamState::SetQuadSolution(QuadSolution type) } break; default: - cout << "Unknon quadrarture function representation" << endl; + cout << "Unknown quadrature function representation" << endl; return; } From 597a94479daebdc8e5aa9a7cc5812c4eec78bde6 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 9 Jul 2024 17:10:51 -0700 Subject: [PATCH 38/57] Added L2 projection. --- CHANGELOG | 9 ++++---- README.md | 3 ++- lib/stream_reader.cpp | 52 +++++++++++++++++++++++++++++++++++++++---- lib/stream_reader.hpp | 3 ++- 4 files changed, 57 insertions(+), 10 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index 968b0dbc..a155d5ba 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -23,11 +23,12 @@ Version 4.2.1 (development) - Added visualization of quadratures. They can be loaded from a file through the new command line argument '-q' or in a socket stream with the keyword 'quadrature' - (instead of solution). Two options of visualization are offered (which can be + (instead of solution). Three options of visualization are offered (which can be switched by 'Q' key): piece-wise constant function on a refined mesh (closed - Gauss-Legendre refinement), or discontinuous elements with DOFs collocated with - the quadrature on the original mesh. High-order quadratures are supported only - for tensor finite elements. + Gauss-Legendre refinement), discontinuous elements with DOFs collocated with + the quadrature on the original mesh, or projection to discontinuous elements. + High-order quadratures are supported only for tensor finite elements with the + first two options. Version 4.2 released on May 23, 2022 diff --git a/README.md b/README.md index 3c6c806f..8e9443ae 100644 --- a/README.md +++ b/README.md @@ -144,7 +144,8 @@ Key commands - Ctrl + p – Print to a PDF file using `gl2ps`. Other vector formats (SVG, EPS) are also possible, but keep in mind that the printing takes a while and the generated files are big. -- Q – Cycle between representations of the quadrature (piece-wise constant refined/L2 element dof collocation) +- Q – Cycle between representations of the quadrature (piece-wise constant refined + / L2 element dof collocation / L2 element projection) - q – Exit ## Advanced diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 8118d2ab..17180a31 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -223,14 +223,18 @@ void StreamState::SetQuadSolution(QuadSolution type) } //check for tensor-product basis - if (order > 0) + if (order > 0 && type != QuadSolution::HO_L2_projected) { Array geoms; mesh->GetGeometries(mesh->Dimension(), geoms); for (auto geom : geoms) if (!Geometry::IsTensorProduct(geom)) { - MFEM_ABORT("High-order quadratures are available only for tensor-product finite elements"); + cout << "High-order quadratures are available only for " + << "tensor-product finite elements with this representation" + << endl; + SetQuadSolution(QuadSolution::HO_L2_projected); + return; } } @@ -241,7 +245,7 @@ void StreamState::SetQuadSolution(QuadSolution type) const int ref_factor = order + 1; if (ref_factor <= 1) { - SetQuadSolution(QuadSolution::HO_L2);//low-order + SetQuadSolution(QuadSolution::HO_L2_collocated);//low-order return; } Mesh *mesh_lor = new Mesh(Mesh::MakeRefined(*mesh, ref_factor, @@ -259,7 +263,7 @@ void StreamState::SetQuadSolution(QuadSolution type) internal.mesh.reset(mesh_lor); } break; - case QuadSolution::HO_L2: + case QuadSolution::HO_L2_collocated: { FiniteElementCollection *fec = new L2_FECollection(order, mesh->Dimension()); FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, @@ -272,6 +276,46 @@ void StreamState::SetQuadSolution(QuadSolution type) internal.grid_f.reset(gf); } break; + case QuadSolution::HO_L2_projected: + { + FiniteElementCollection *fec = new L2_FECollection(order, mesh->Dimension()); + FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, + quad_f->GetVDim(), Ordering::byVDIM); + cout << "Projecting quadrature to grid function of order " + << order << endl; + GridFunction *gf = new GridFunction(fes); + gf->MakeOwner(fec); + + VectorQuadratureFunctionCoefficient coeff_qf(*quad_f); + VectorDomainLFIntegrator bi(coeff_qf); + VectorMassIntegrator Mi; + Mi.SetVDim(quad_f->GetVDim()); + + DenseMatrix M_i; + Vector x_i, b_i; + DenseMatrixInverse M_i_inv; + Array vdofs; + + for (int i = 0; i < mesh->GetNE(); i++) + { + const FiniteElement *fe = fes->GetFE(i); + ElementTransformation *Tr = mesh->GetElementTransformation(i); + Mi.AssembleElementMatrix(*fe, *Tr, M_i); + M_i_inv.Factor(M_i); + + bi.SetIntRule(&quad_f->GetIntRule(i)); + bi.AssembleRHSElementVect(*fe, *Tr, b_i); + + x_i.SetSize(b_i.Size()); + M_i_inv.Mult(b_i, x_i); + + fes->GetElementVDofs(i, vdofs); + gf->SetSubVector(vdofs, x_i); + } + + internal.grid_f.reset(gf); + } + break; default: cout << "Unknown quadrature function representation" << endl; return; diff --git a/lib/stream_reader.hpp b/lib/stream_reader.hpp index ab02f73c..4c004aa9 100644 --- a/lib/stream_reader.hpp +++ b/lib/stream_reader.hpp @@ -63,7 +63,8 @@ struct StreamState MIN = -1, //---------- LOR_ClosedGL, - HO_L2, + HO_L2_collocated, + HO_L2_projected, //---------- MAX } quad_sol {QuadSolution::NONE}; From e026126e8d26514ad63b89693e71f92ce1d0162f Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 17 Jul 2024 11:30:13 -0700 Subject: [PATCH 39/57] Changed the word quadrature in CHANGELOG --- CHANGELOG | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index a155d5ba..4d217d7b 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -21,14 +21,14 @@ Version 4.2.1 (development) - Added a compilation parameter for the default font size. -- Added visualization of quadratures. They can be loaded from a file through the - new command line argument '-q' or in a socket stream with the keyword 'quadrature' - (instead of solution). Three options of visualization are offered (which can be - switched by 'Q' key): piece-wise constant function on a refined mesh (closed - Gauss-Legendre refinement), discontinuous elements with DOFs collocated with - the quadrature on the original mesh, or projection to discontinuous elements. - High-order quadratures are supported only for tensor finite elements with the - first two options. +- Added visualization of quadrature data (QuadratureFunction in MFEM). They can + be loaded from a file through the new command line argument '-q' or in a socket + stream with the keyword 'quadrature' (instead of 'solution'). Three options of + visualization are offered (which can be switched by 'Q' key): piece-wise constant + function on a refined mesh (closed Gauss-Legendre refinement), discontinuous + elements with DOFs collocated with the quadrature on the original mesh, or + projection to discontinuous elements. High-order quadrature functions are supported + only for tensor finite elements with the first two options. Version 4.2 released on May 23, 2022 From 226fd1bafe50a2d8b50a1ef8a93a76be607cc8ab Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 17 Jul 2024 11:33:42 -0700 Subject: [PATCH 40/57] Improved README.md --- README.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 8e9443ae..3795d329 100644 --- a/README.md +++ b/README.md @@ -144,8 +144,10 @@ Key commands - Ctrl + p – Print to a PDF file using `gl2ps`. Other vector formats (SVG, EPS) are also possible, but keep in mind that the printing takes a while and the generated files are big. -- Q – Cycle between representations of the quadrature (piece-wise constant refined - / L2 element dof collocation / L2 element projection) +- Q – Cycle between representations of the visualized *quadrature data*. The options are: + - piece-wise constant refined + - L2 element dof collocation + - L2 element projection - q – Exit ## Advanced From 137c9a3b1f7bf23e3422acd7610a5443017d021d Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 17 Jul 2024 11:58:36 -0700 Subject: [PATCH 41/57] Added Q key mention to the help strings. --- lib/vssolution.cpp | 1 + lib/vssolution3d.cpp | 1 + lib/vsvector.cpp | 1 + lib/vsvector3d.cpp | 1 + 4 files changed, 4 insertions(+) diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 0fc416a4..a6861b59 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -81,6 +81,7 @@ std::string VisualizationSceneSolution::GetHelpString() const << "| s - Turn on/off unit cube scaling |" << endl << "| S - Take snapshot/Record a movie |" << endl << "| t - Cycle materials and lights |" << endl + << "| Q - Cycle quadrature data repre. |" << endl << "| y/Y Rotate the cutting plane |" << endl << "| z/Z Move the cutting plane |" << endl << "| \\ - Set light source position |" << endl diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 47c543ec..7ed7de97 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -69,6 +69,7 @@ std::string VisualizationSceneSolution3d::GetHelpString() const << "| s - Turn on/off unit cube scaling |" << endl << "| S - Take snapshot/Record a movie |" << endl << "| t - Cycle materials and lights |" << endl + << "| Q - Cycle quadrature data repre. |" << endl << "| u/U Move the level surface |" << endl << "| v/V Add/Delete a level surface |" << endl << "| w/W Move bdr elements up/down |" << endl diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index db133da9..c83b8aa1 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -56,6 +56,7 @@ std::string VisualizationSceneVector::GetHelpString() const << "| s - Turn on/off unit cube scaling |" << endl << "| S - Take snapshot/Record a movie |" << endl << "| t - Cycle materials and lights |" << endl + << "| Q - Cycle quadrature data repre. |" << endl << "| u - Vector sampling; scalar func. |" << endl << "| U - Switch 'u' functionality |" << endl << "| v - Cycle through vector fields |" << endl diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index 2a41b299..a081c433 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -66,6 +66,7 @@ std::string VisualizationSceneVector3d::GetHelpString() const << "| s - Turn on/off unit cube scaling |" << endl << "| S - Take snapshot/Record a movie |" << endl << "| t - Cycle materials and lights |" << endl + << "| Q - Cycle quadrature data repre. |" << endl << "| u/U Move the level field vectors |" << endl << "| v/V Vector field |" << endl << "| w/W Add/Delete level field vector |" << endl From 6b6c6f3736adb4fec76020a888eb75ef7b53dbb8 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 17 Jul 2024 12:01:20 -0700 Subject: [PATCH 42/57] Changed quadrature -> quadrature data. --- lib/stream_reader.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 17180a31..4c8e657c 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -230,7 +230,7 @@ void StreamState::SetQuadSolution(QuadSolution type) for (auto geom : geoms) if (!Geometry::IsTensorProduct(geom)) { - cout << "High-order quadratures are available only for " + cout << "High-order quadrature data are supported only for " << "tensor-product finite elements with this representation" << endl; SetQuadSolution(QuadSolution::HO_L2_projected); @@ -254,7 +254,7 @@ void StreamState::SetQuadSolution(QuadSolution type) FiniteElementSpace *fes = new FiniteElementSpace(mesh_lor, fec, quad_f->GetVDim(), Ordering::byVDIM); MFEM_ASSERT(quad_f->Size() == fes->GetVSize(), "Size mismatch"); - cout << "Representing quadrature by piecewise-constant function on mesh refined " + cout << "Representing quadrature data by piecewise-constant function on mesh refined " << ref_factor << " times" << endl; GridFunction *gf = new GridFunction(fes, *quad_f, 0); gf->MakeOwner(fec); @@ -269,7 +269,7 @@ void StreamState::SetQuadSolution(QuadSolution type) FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, quad_f->GetVDim(), Ordering::byVDIM); MFEM_ASSERT(quad_f->Size() == fes->GetVSize(), "Size mismatch"); - cout << "Representing quadrature by grid function of order " + cout << "Representing quadrature data by grid function of order " << order << endl; GridFunction *gf = new GridFunction(fes, *quad_f, 0); gf->MakeOwner(fec); @@ -281,7 +281,7 @@ void StreamState::SetQuadSolution(QuadSolution type) FiniteElementCollection *fec = new L2_FECollection(order, mesh->Dimension()); FiniteElementSpace *fes = new FiniteElementSpace(mesh.get(), fec, quad_f->GetVDim(), Ordering::byVDIM); - cout << "Projecting quadrature to grid function of order " + cout << "Projecting quadrature data to grid function of order " << order << endl; GridFunction *gf = new GridFunction(fes); gf->MakeOwner(fec); @@ -317,7 +317,7 @@ void StreamState::SetQuadSolution(QuadSolution type) } break; default: - cout << "Unknown quadrature function representation" << endl; + cout << "Unknown quadrature data representation" << endl; return; } From 47c1935cbacf83e415a860e2f429a1ddaeab53a1 Mon Sep 17 00:00:00 2001 From: Tzanio Kolev Date: Wed, 17 Jul 2024 12:24:38 -0700 Subject: [PATCH 43/57] minor --- CHANGELOG | 15 +++++++-------- README.md | 6 +++--- lib/vssolution.cpp | 2 +- lib/vssolution3d.cpp | 2 +- lib/vsvector.cpp | 2 +- lib/vsvector3d.cpp | 2 +- 6 files changed, 14 insertions(+), 15 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index d1cd43b8..8010c9d7 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -25,14 +25,13 @@ Version 4.2.1 (development) - Added a compilation parameter for the default font size. -- Added visualization of quadrature data (QuadratureFunction in MFEM). They can - be loaded from a file through the new command line argument '-q' or in a socket - stream with the keyword 'quadrature' (instead of 'solution'). Three options of - visualization are offered (which can be switched by 'Q' key): piece-wise constant - function on a refined mesh (closed Gauss-Legendre refinement), discontinuous - elements with DOFs collocated with the quadrature on the original mesh, or - projection to discontinuous elements. High-order quadrature functions are supported - only for tensor finite elements with the first two options. +- Added visualization of quadrature data (QuadratureFunction in MFEM). Both + loading from file, with the new command line argument '-q', or from a socket + stream, with the keyword 'quadrature', are supported. Three visualization + options are provided: piece-wise constants on a refined mesh (LOR), L2 field + with DOFs collocated (interpolation), or projection to discontinuous elements + (L2 projection). Use 'Q' to switch between them. High-order quadrature data is + supported only for tensor finite elements with the first two options. Version 4.2 released on May 23, 2022 diff --git a/README.md b/README.md index 3795d329..ca73cc36 100644 --- a/README.md +++ b/README.md @@ -145,9 +145,9 @@ Key commands vector formats (SVG, EPS) are also possible, but keep in mind that the printing takes a while and the generated files are big. - Q – Cycle between representations of the visualized *quadrature data*. The options are: - - piece-wise constant refined - - L2 element dof collocation - - L2 element projection + - piece-wise constant refined (LOR) + - L2 element dof collocation (interpolation) + - L2 element projection (L2 projection) - q – Exit ## Advanced diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 2eacaf97..b29e2de2 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -76,12 +76,12 @@ std::string VisualizationSceneSolution::GetHelpString() const << "| O - Switch 'o' func. (NC shading) |" << endl << "| p/P Cycle through color palettes |" << endl << "| q - Quits |" << endl + << "| Q - Cycle quadrature data mode |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl << "| S - Take snapshot/Record a movie |" << endl << "| t - Cycle materials and lights |" << endl - << "| Q - Cycle quadrature data repre. |" << endl << "| y/Y Rotate the cutting plane |" << endl << "| z/Z Move the cutting plane |" << endl << "| \\ - Set light source position |" << endl diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 4493bd48..41d6cb76 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -64,12 +64,12 @@ std::string VisualizationSceneSolution3d::GetHelpString() const << "| o/O (De)refine elem, disc shading |" << endl << "| p/P Cycle through color palettes |" << endl << "| q - Quits |" << endl + << "| Q - Cycle quadrature data mode |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl << "| S - Take snapshot/Record a movie |" << endl << "| t - Cycle materials and lights |" << endl - << "| Q - Cycle quadrature data repre. |" << endl << "| u/U Move the level surface |" << endl << "| v/V Add/Delete a level surface |" << endl << "| w/W Move bdr elements up/down |" << endl diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index c83b8aa1..909f5396 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -51,12 +51,12 @@ std::string VisualizationSceneVector::GetHelpString() const << "| O - Switch 'o' func. (NC shading) |" << endl << "| p/P Cycle through color palettes |" << endl << "| q - Quits |" << endl + << "| Q - Cycle quadrature data mode |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl << "| S - Take snapshot/Record a movie |" << endl << "| t - Cycle materials and lights |" << endl - << "| Q - Cycle quadrature data repre. |" << endl << "| u - Vector sampling; scalar func. |" << endl << "| U - Switch 'u' functionality |" << endl << "| v - Cycle through vector fields |" << endl diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index a081c433..4f5759d3 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -61,12 +61,12 @@ std::string VisualizationSceneVector3d::GetHelpString() const << "| o/O (De)refine elem, disc shading |" << endl << "| p/P Cycle through color palettes |" << endl << "| q - Quits |" << endl + << "| Q - Cycle quadrature data mode |" << endl << "| r - Reset the plot to 3D view |" << endl << "| R - Reset the plot to 2D view |" << endl << "| s - Turn on/off unit cube scaling |" << endl << "| S - Take snapshot/Record a movie |" << endl << "| t - Cycle materials and lights |" << endl - << "| Q - Cycle quadrature data repre. |" << endl << "| u/U Move the level field vectors |" << endl << "| v/V Vector field |" << endl << "| w/W Add/Delete level field vector |" << endl From f25e99eb03776702991b9d88ac95cb44a1282c61 Mon Sep 17 00:00:00 2001 From: Tzanio Kolev Date: Thu, 18 Jul 2024 09:59:11 -0700 Subject: [PATCH 44/57] updated test/data --- tests/data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/data b/tests/data index 53fbee22..81de802a 160000 --- a/tests/data +++ b/tests/data @@ -1 +1 @@ -Subproject commit 53fbee222cceb0581122cb72518bfc1b85766b4e +Subproject commit 81de802ad63ba41cd1fcd3825fbb3d2471538ace From 8c2dad09e9ed5707d827e57157b46efb9c455265 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 10:42:50 -0700 Subject: [PATCH 45/57] Changed shading from int to enum class. --- glvis.cpp | 13 ++++---- lib/threads.cpp | 11 ++++--- lib/vsdata.hpp | 14 ++++++++- lib/vssolution.cpp | 51 +++++++++++++++--------------- lib/vssolution.hpp | 6 ++-- lib/vssolution3d.cpp | 74 +++++++++++++++++++++++++------------------- lib/vssolution3d.hpp | 7 +++-- lib/vsvector.cpp | 9 +++--- lib/vsvector3d.cpp | 6 ++-- 9 files changed, 109 insertions(+), 82 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index f9d41508..e48a535f 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -255,7 +255,7 @@ bool GLVisInitVis(StreamState::FieldType field_type, if (stream_state.grid_f) { vs->AutoRefine(); - vs->SetShading(2, true); + vs->SetShading(VisualizationSceneScalarData::Shading::Noncomforming, true); } if (mesh_range > 0.0) { @@ -686,20 +686,21 @@ void ExecuteScriptCommand() { scr >> ws >> word; cout << "Script: shading: " << flush; - int s = -1; + VisualizationSceneScalarData::Shading s = + VisualizationSceneScalarData::Shading::Invalid; if (word == "flat") { - s = 0; + s = VisualizationSceneScalarData::Shading::Flat; } else if (word == "smooth") { - s = 1; + s = VisualizationSceneScalarData::Shading::Smooth; } else if (word == "cool") { - s = 2; + s = VisualizationSceneScalarData::Shading::Noncomforming; } - if (s != -1) + if (s != VisualizationSceneScalarData::Shading::Invalid) { vs->SetShading(s, false); cout << word << endl; diff --git a/lib/threads.cpp b/lib/threads.cpp index 23d309fe..9044fc19 100644 --- a/lib/threads.cpp +++ b/lib/threads.cpp @@ -580,20 +580,21 @@ int GLVisCommand::Execute() case SHADING: { cout << "Command: shading: " << flush; - int s = -1; + VisualizationSceneScalarData::Shading s = + VisualizationSceneScalarData::Shading::Invalid; if (shading == "flat") { - s = 0; + s = VisualizationSceneScalarData::Shading::Flat; } else if (shading == "smooth") { - s = 1; + s = VisualizationSceneScalarData::Shading::Smooth; } else if (shading == "cool") { - s = 2; + s = VisualizationSceneScalarData::Shading::Noncomforming; } - if (s != -1) + if (s != VisualizationSceneScalarData::Shading::Invalid) { (*vs)->SetShading(s, false); cout << shading << endl; diff --git a/lib/vsdata.hpp b/lib/vsdata.hpp index d08fedd6..8214cea8 100644 --- a/lib/vsdata.hpp +++ b/lib/vsdata.hpp @@ -127,6 +127,18 @@ class VisualizationSceneScalarData : public VisualizationScene void Cone(gl3::GlDrawable& buf, glm::mat4 transform, double cval); public: + enum class Shading + { + Invalid = -1, + Min = -1, + //--------- + Flat, + Smooth, + Noncomforming, + //--------- + Max + }; + Plane *CuttingPlane; int key_r_state; /** Shrink factor with respect to the center of each element (2D) or the @@ -182,7 +194,7 @@ class VisualizationSceneScalarData : public VisualizationScene virtual void UpdateValueRange(bool prepare) = 0; void SetValueRange(double, double); - virtual void SetShading(int, bool) = 0; + virtual void SetShading(Shading, bool) = 0; virtual void SetRefineFactors(int, int) = 0; void SetAutoRefineLimits(int max_ref, int max_surf_elem) { diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index b29e2de2..c81ad650 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -356,7 +356,7 @@ static void KeyZPressed() static void KeyF3Pressed() { - if (vssol->shading == 2) + if (vssol->shading == VisualizationSceneScalarData::Shading::Noncomforming) { vssol->shrink *= 0.9; vssol->Prepare(); @@ -370,7 +370,7 @@ static void KeyF3Pressed() static void KeyF4Pressed() { - if (vssol->shading == 2) + if (vssol->shading == VisualizationSceneScalarData::Shading::Noncomforming) { vssol->shrink *= 1.11111111111111111111111; vssol->Prepare(); @@ -383,7 +383,7 @@ static void KeyF4Pressed() static void KeyF11Pressed() { - if (vssol->shading == 2) + if (vssol->shading == VisualizationSceneScalarData::Shading::Noncomforming) { if (vssol->matc.Width() == 0) { @@ -401,7 +401,7 @@ static void KeyF11Pressed() static void KeyF12Pressed() { - if (vssol->shading == 2) + if (vssol->shading == VisualizationSceneScalarData::Shading::Noncomforming) { if (vssol->matc.Width() == 0) { @@ -437,7 +437,8 @@ void VisualizationSceneSolution::Init() rsol = NULL; vssol = this; - drawelems = shading = 1; + drawelems = 1; + shading = Shading::Smooth; drawmesh = 0; draworder = 0; drawnums = 0; @@ -558,7 +559,7 @@ void VisualizationSceneSolution::ToggleDrawElems() maxv_sol = maxv; have_sol_range = true; } - else if (shading == 2) + else if (shading == Shading::Noncomforming) { if (drawelems == 1 && have_sol_range) { @@ -823,21 +824,21 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals( return have_normals; } -void VisualizationSceneSolution::SetShading(int s, bool print) +void VisualizationSceneSolution::SetShading(Shading s, bool print) { - if (shading == s || s < 0) + if (shading == s || s <= Shading::Min) { return; } if (rsol) { - if (s > 2) + if (s >= Shading::Max) { return; } - if (s == 2 || shading == 2) + if (s == Shading::Noncomforming || shading == Shading::Noncomforming) { shading = s; have_sol_range = false; @@ -856,7 +857,7 @@ void VisualizationSceneSolution::SetShading(int s, bool print) } else { - if (s > 1) + if (s > Shading::Smooth) { return; } @@ -868,7 +869,7 @@ void VisualizationSceneSolution::SetShading(int s, bool print) {"flat", "smooth", "non-conforming (with subdivision)"}; if (print) { - cout << "Shading type : " << shading_type[shading] << endl; + cout << "Shading type : " << shading_type[(int)shading] << endl; } } @@ -876,11 +877,11 @@ void VisualizationSceneSolution::ToggleShading() { if (rsol) { - SetShading((shading + 1) % 3, true); + SetShading((Shading)(((int)shading + 1) % (int)Shading::Max), true); } else { - SetShading(1 - shading, true); + SetShading((Shading)(1 - (int)shading), true); } } @@ -920,7 +921,7 @@ void VisualizationSceneSolution::ToggleRefinements() } break; } - if (update && shading == 2) + if (update && shading == Shading::Noncomforming) { have_sol_range = false; DoAutoscale(false); @@ -971,7 +972,7 @@ void VisualizationSceneSolution::SetRefineFactors(int tot, int bdr) TimesToRefine = tot; EdgeRefineFactor = bdr; - if (shading == 2) + if (shading == Shading::Noncomforming) { have_sol_range = false; DoAutoscale(false); @@ -1054,7 +1055,7 @@ void VisualizationSceneSolution::FindNewBox(double rx[], double ry[], { int i, j; - if (shading != 2) + if (shading != Shading::Noncomforming) { int nv = mesh -> GetNV(); @@ -1437,10 +1438,10 @@ void VisualizationSceneSolution::Prepare() switch (shading) { - case 0: + case Shading::Flat: PrepareFlat(); return; - case 2: + case Shading::Noncomforming: PrepareFlat2(); return; default: @@ -1553,7 +1554,7 @@ void VisualizationSceneSolution::Prepare() void VisualizationSceneSolution::PrepareLevelCurves() { - if (shading == 2) + if (shading == Shading::Noncomforming) { PrepareLevelCurves2(); return; @@ -1691,7 +1692,7 @@ void VisualizationSceneSolution::PrepareLevelCurves2() void VisualizationSceneSolution::PrepareLines() { - if (shading == 2 && + if (shading == Shading::Noncomforming && mesh->Dimension() > 1) // PrepareLines3 does not make sense for 1d meshes. { // PrepareLines2(); @@ -1762,7 +1763,7 @@ double VisualizationSceneSolution::GetElementLengthScale(int k) void VisualizationSceneSolution::PrepareElementNumbering() { - if (2 == shading) + if (shading == Shading::Noncomforming) { PrepareElementNumbering2(); } @@ -1844,7 +1845,7 @@ void VisualizationSceneSolution::PrepareElementNumbering2() void VisualizationSceneSolution::PrepareVertexNumbering() { - if (2 == shading) + if (shading == Shading::Noncomforming) { PrepareVertexNumbering2(); } @@ -2190,7 +2191,7 @@ void VisualizationSceneSolution::PrepareBoundary() bdr_buf.clear(); gl3::GlBuilder bl = bdr_buf.createBuilder(); - if (shading != 2) + if (shading != Shading::Noncomforming) { bl.glBegin(GL_LINES); for (i = 0; i < ne; i++) @@ -2277,7 +2278,7 @@ void VisualizationSceneSolution::PrepareCP() gl3::GlBuilder bld = cp_buf.createBuilder(); bld.glBegin(GL_LINES); - if (shading != 2) + if (shading != Shading::Noncomforming) { Array vertices; diff --git a/lib/vssolution.hpp b/lib/vssolution.hpp index c86cfc2c..b1cd5e33 100644 --- a/lib/vssolution.hpp +++ b/lib/vssolution.hpp @@ -83,7 +83,7 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData double GetElementLengthScale(int k); public: - int shading; + Shading shading; int attr_to_show, bdr_attr_to_show; Array el_attr_to_show, bdr_el_attr_to_show; @@ -165,7 +165,7 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData PrepareNumbering(false); } - virtual void SetShading(int, bool); + virtual void SetShading(Shading, bool); void ToggleShading(); void ToggleDrawCP() { draw_cp = !draw_cp; PrepareCP(); } @@ -178,7 +178,7 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData virtual void ToggleAttributes(Array &attr_list); virtual void SetDrawMesh(int i) { drawmesh = i % 3; } - virtual int GetShading() { return shading; } + virtual Shading GetShading() { return shading; } virtual int GetDrawMesh() { return drawmesh; } }; diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 41d6cb76..508eddf8 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -349,7 +349,8 @@ static void KeyoPressed(GLenum state) if (vssol3d -> TimesToRefine < 32) { cout << "Subdivision factor = " << ++vssol3d->TimesToRefine << endl; - if (vssol3d -> GetShading() == 2) + if (vssol3d -> GetShading() == + VisualizationSceneScalarData::Shading::Noncomforming) { vssol3d->DoAutoscale(false); vssol3d -> Prepare(); @@ -367,7 +368,8 @@ static void KeyOPressed() if (vssol3d -> TimesToRefine > 1) { cout << "Subdivision factor = " << --vssol3d->TimesToRefine << endl; - if (vssol3d -> GetShading() == 2) + if (vssol3d -> GetShading() == + VisualizationSceneScalarData::Shading::Noncomforming) { vssol3d->DoAutoscale(false); vssol3d -> Prepare(); @@ -381,7 +383,8 @@ static void KeyOPressed() static void KeywPressed() { - if (vssol3d -> GetShading() == 2) + if (vssol3d -> GetShading() == + VisualizationSceneScalarData::Shading::Noncomforming) { if ( fabs(vssol3d -> FaceShiftScale += 0.01) < 0.001 ) { @@ -398,7 +401,8 @@ static void KeywPressed() static void KeyWPressed() { - if (vssol3d -> GetShading() == 2) + if (vssol3d -> GetShading() == + VisualizationSceneScalarData::Shading::Noncomforming) { if ( fabs(vssol3d -> FaceShiftScale -= 0.01) < 0.001 ) { @@ -461,7 +465,8 @@ static void KeyF3Pressed(GLenum state) } else { - if (vssol3d->GetShading() == 2) + if (vssol3d->GetShading() == + VisualizationSceneScalarData::Shading::Noncomforming) { if (vssol3d->GetMesh()->Dimension() == 3 && vssol3d->bdrc.Width() == 0) { @@ -502,7 +507,8 @@ static void KeyF4Pressed(GLenum state) } else { - if (vssol3d->GetShading() == 2) + if (vssol3d->GetShading() == + VisualizationSceneScalarData::Shading::Noncomforming) { if (vssol3d->GetMesh()->Dimension() == 3 && vssol3d->bdrc.Width() == 0) { @@ -526,7 +532,8 @@ static void KeyF4Pressed(GLenum state) static void KeyF11Pressed() { - if (vssol3d->GetShading() == 2) + if (vssol3d->GetShading() == + VisualizationSceneScalarData::Shading::Noncomforming) { if (vssol3d->matc.Width() == 0) { @@ -545,7 +552,8 @@ static void KeyF11Pressed() static void KeyF12Pressed() { - if (vssol3d->GetShading() == 2) + if (vssol3d->GetShading() == + VisualizationSceneScalarData::Shading::Noncomforming) { if (vssol3d->matc.Width() == 0) { @@ -699,7 +707,8 @@ void VisualizationSceneSolution3d::Init() drawlsurf = 0; cp_algo = 0; - drawelems = shading = 1; + drawelems = 1; + shading = Shading::Smooth; drawmesh = 0; draworder = 0; scaling = 0; @@ -843,20 +852,21 @@ void VisualizationSceneSolution3d::NewMeshAndSolution( PrepareOrderingCurve(); } -void VisualizationSceneSolution3d::SetShading(int s, bool print) +void VisualizationSceneSolution3d::SetShading(Shading s, bool print) { - if (shading == s || s < 0) + if (shading == s || s <= Shading::Min) { return; } - if (s > 2 || (GridF == NULL && s > 1)) + if (s >= Shading::Max || (GridF == NULL && s > Shading::Smooth)) { return; } - int os = shading; + Shading os = shading; shading = s; - if (GridF != NULL && (s == 2 || os == 2)) + if (GridF != NULL && (s == Shading::Noncomforming || + os == Shading::Noncomforming)) { DoAutoscale(false); PrepareLines(); @@ -869,7 +879,7 @@ void VisualizationSceneSolution3d::SetShading(int s, bool print) {"flat", "smooth", "non-conforming (with subdivision)"}; if (print) { - cout << "Shading type : " << shading_type[shading] << endl; + cout << "Shading type : " << shading_type[(int)shading] << endl; } } @@ -877,11 +887,11 @@ void VisualizationSceneSolution3d::ToggleShading() { if (GridF) { - SetShading((shading+1)%3, true); + SetShading((Shading)(((int)shading+1) % (int)Shading::Max), true); } else { - SetShading(1-shading, true); + SetShading((Shading)(1 - (int)shading), true); } } @@ -894,7 +904,7 @@ void VisualizationSceneSolution3d::SetRefineFactors(int f, int ignored) TimesToRefine = f; - if (shading == 2) + if (shading == Shading::Noncomforming) { DoAutoscale(false); Prepare(); @@ -979,7 +989,7 @@ void VisualizationSceneSolution3d::FindNewBox(bool prepare) if (coord[2] > bb.z[1]) { bb.z[1] = coord[2]; } } - if (shading == 2) + if (shading == Shading::Noncomforming) { int dim = mesh->Dimension(); int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE(); @@ -1030,7 +1040,7 @@ void VisualizationSceneSolution3d::FindNewValueRange(bool prepare) GridF->FESpace()->FEColl()->GetMapType(mesh->Dimension()) : FiniteElement::VALUE; - if (shading < 2 || map_type != (int)FiniteElement::VALUE) + if (shading < Shading::Noncomforming || map_type != (int)FiniteElement::VALUE) { minv = sol->Min(); maxv = sol->Max(); @@ -1050,7 +1060,7 @@ void VisualizationSceneSolution3d::EventUpdateColors() PrepareCuttingPlane(); PrepareLevelSurf(); PrepareOrderingCurve(); - if (shading == 2 && drawmesh != 0 && FaceShiftScale != 0.0) + if (shading == Shading::Noncomforming && drawmesh != 0 && FaceShiftScale != 0.0) { PrepareLines(); } @@ -1127,7 +1137,7 @@ void VisualizationSceneSolution3d::ToggleCPDrawMesh() void VisualizationSceneSolution3d::ToggleCPAlgorithm() { cp_algo = (cp_algo+1)%2; - if (shading == 2 && cplane == 1) + if (shading == Shading::Noncomforming && cplane == 1) { CPPrepare(); } @@ -1817,10 +1827,10 @@ void VisualizationSceneSolution3d::Prepare() switch (shading) { - case 0: + case Shading::Flat: PrepareFlat(); return; - case 2: + case Shading::Noncomforming: PrepareFlat2(); return; default: @@ -2003,7 +2013,7 @@ void VisualizationSceneSolution3d::PrepareLines() return; } - if (shading == 2) + if (shading == Shading::Noncomforming) { PrepareLines2(); return; @@ -2523,7 +2533,7 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc(int func) while (n > 2) { - if (shading != 2) + if (shading != Shading::Noncomforming) { mesh -> GetPointMatrix (i, pointmat); } @@ -2556,7 +2566,7 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc(int func) { case 1: // PrepareCuttingPlane() { - if (shading == 2) + if (shading == Shading::Noncomforming) { // changes point for n > 4 DrawRefinedSurf(n, point[0], i, 1); @@ -2606,7 +2616,7 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc(int func) case 2: // PrepareCuttingPlaneLines() with mesh { - if (shading == 2) + if (shading == Shading::Noncomforming) { // changes point for n > 4 DrawRefinedSurf(n, point[0], i, 2); @@ -2627,7 +2637,7 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc(int func) case 3: // PrepareCuttingPlaneLines() with level lines { - if (shading == 2) + if (shading == Shading::Noncomforming) { // changes point for n > 4 DrawRefinedSurf(n, point[0], i, 3); @@ -2918,7 +2928,7 @@ void VisualizationSceneSolution3d::PrepareCuttingPlane2() mesh -> GetFaceElements (i, &e1, &e2); if (e2 >= 0 && partition[e1] != partition[e2]) { - if (shading != 2) + if (shading != Shading::Noncomforming) { mesh -> GetFaceVertices (i, nodes); for (j = 0; j < nodes.Size(); j++) @@ -3080,7 +3090,7 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() mesh -> GetFaceElements (i, &e1, &e2); if (e2 >= 0 && partition[e1] != partition[e2]) { - if (shading != 2) + if (shading != Shading::Noncomforming) { mesh -> GetFaceVertices (i, nodes); for (j = 0; j < nodes.Size(); j++) @@ -3869,7 +3879,7 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() Array faces, ofaces; - if (shading != 2) + if (shading != Shading::Noncomforming) { for (int ie = 0; ie < mesh->GetNE(); ie++) { diff --git a/lib/vssolution3d.hpp b/lib/vssolution3d.hpp index a184e403..2179baf2 100644 --- a/lib/vssolution3d.hpp +++ b/lib/vssolution3d.hpp @@ -22,7 +22,8 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData { protected: - int drawmesh, drawelems, shading, draworder; + int drawmesh, drawelems, draworder; + Shading shading; int cplane; int cp_drawmesh, cp_drawelems, drawlsurf; // Algorithm used to draw the cutting plane when shading is 2 and cplane is 1 @@ -148,8 +149,8 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData void ToggleDrawOrdering() { draworder = (draworder+1)%5; } void ToggleShading(); - int GetShading() { return shading; }; - virtual void SetShading(int, bool); + Shading GetShading() { return shading; }; + virtual void SetShading(Shading, bool); virtual void SetRefineFactors(int, int); virtual void AutoRefine(); virtual void ToggleAttributes(Array &attr_list); diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index 909f5396..0149b6c1 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -187,7 +187,8 @@ void KeyuPressed() { case 0: case 1: - if (update && vsvector->shading == 2) + if (update && + vsvector->shading == VisualizationSceneSolution::Shading::Noncomforming) { vsvector->PrepareVectorField(); SendExposeEvent(); @@ -230,7 +231,7 @@ void VisualizationSceneVector::ToggleDrawElems() cout << "Surface elements mode : " << modes[drawelems] << endl; - if (drawelems != 0 && shading == 2) + if (drawelems != 0 && shading == Shading::Noncomforming) { DoAutoscaleValue(false); PrepareLines(); @@ -660,7 +661,7 @@ void VisualizationSceneVector::PrepareDisplacedMesh() // prepare the displaced mesh displine_buf.clear(); gl3::GlBuilder build = displine_buf.createBuilder(); - if (shading != 2) + if (shading != Shading::Noncomforming) { for (int i = 0; i < ne; i++) { @@ -905,7 +906,7 @@ void VisualizationSceneVector::PrepareVectorField() DrawVector(v[0], v[1], (*solx)(i), (*soly)(i), (*sol)(i)); } - if (shading == 2 && RefineFactor > 1) + if (shading == Shading::Noncomforming && RefineFactor > 1) { DenseMatrix vvals, pm; for (i = 0; i < mesh->GetNE(); i++) diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index 4f5759d3..623a7497 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -794,10 +794,10 @@ void VisualizationSceneVector3d::Prepare() switch (shading) { - case 0: + case Shading::Flat: PrepareFlat(); return; - case 2: + case Shading::Noncomforming: PrepareFlat2(); return; default: @@ -928,7 +928,7 @@ void VisualizationSceneVector3d::PrepareLines() { if (!drawmesh) { return; } - if (shading == 2) + if (shading == Shading::Noncomforming) { PrepareLines2(); return; From c730066f37d398241cac5b63affa810821be712e Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 10:56:17 -0700 Subject: [PATCH 46/57] Moved shading to VisualizationSceneScalarData. --- lib/vsdata.hpp | 28 ++++++++++++++++------------ lib/vssolution.cpp | 10 +++++----- lib/vssolution.hpp | 5 +---- lib/vssolution3d.cpp | 2 +- lib/vssolution3d.hpp | 3 +-- lib/vsvector.cpp | 2 +- 6 files changed, 25 insertions(+), 25 deletions(-) diff --git a/lib/vsdata.hpp b/lib/vsdata.hpp index 8214cea8..d802f48f 100644 --- a/lib/vsdata.hpp +++ b/lib/vsdata.hpp @@ -54,6 +54,19 @@ class Plane class VisualizationSceneScalarData : public VisualizationScene { +public: + enum class Shading + { + Invalid = -1, + Min = -1, + //--------- + Flat, + Smooth, + Noncomforming, + //--------- + Max + }; + protected: Mesh *mesh; Vector *sol; @@ -63,6 +76,7 @@ class VisualizationSceneScalarData : public VisualizationScene std::string a_label_x, a_label_y, a_label_z; int scaling, colorbar, drawaxes; + Shading shading; int auto_ref_max, auto_ref_max_surf_elem; vector updated_bufs; @@ -127,18 +141,6 @@ class VisualizationSceneScalarData : public VisualizationScene void Cone(gl3::GlDrawable& buf, glm::mat4 transform, double cval); public: - enum class Shading - { - Invalid = -1, - Min = -1, - //--------- - Flat, - Smooth, - Noncomforming, - //--------- - Max - }; - Plane *CuttingPlane; int key_r_state; /** Shrink factor with respect to the center of each element (2D) or the @@ -195,6 +197,8 @@ class VisualizationSceneScalarData : public VisualizationScene void SetValueRange(double, double); virtual void SetShading(Shading, bool) = 0; + virtual void ToogleShading() { SetShading((Shading)(((int)shading + 1) % (int)Shading::Max), true); } + virtual Shading GetShading() { return shading; } virtual void SetRefineFactors(int, int) = 0; void SetAutoRefineLimits(int max_ref, int max_surf_elem) { diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index c81ad650..a3a99270 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -356,7 +356,7 @@ static void KeyZPressed() static void KeyF3Pressed() { - if (vssol->shading == VisualizationSceneScalarData::Shading::Noncomforming) + if (vssol->GetShading() == VisualizationSceneScalarData::Shading::Noncomforming) { vssol->shrink *= 0.9; vssol->Prepare(); @@ -370,7 +370,7 @@ static void KeyF3Pressed() static void KeyF4Pressed() { - if (vssol->shading == VisualizationSceneScalarData::Shading::Noncomforming) + if (vssol->GetShading() == VisualizationSceneScalarData::Shading::Noncomforming) { vssol->shrink *= 1.11111111111111111111111; vssol->Prepare(); @@ -383,7 +383,7 @@ static void KeyF4Pressed() static void KeyF11Pressed() { - if (vssol->shading == VisualizationSceneScalarData::Shading::Noncomforming) + if (vssol->GetShading() == VisualizationSceneScalarData::Shading::Noncomforming) { if (vssol->matc.Width() == 0) { @@ -401,7 +401,7 @@ static void KeyF11Pressed() static void KeyF12Pressed() { - if (vssol->shading == VisualizationSceneScalarData::Shading::Noncomforming) + if (vssol->GetShading() == VisualizationSceneScalarData::Shading::Noncomforming) { if (vssol->matc.Width() == 0) { @@ -877,7 +877,7 @@ void VisualizationSceneSolution::ToggleShading() { if (rsol) { - SetShading((Shading)(((int)shading + 1) % (int)Shading::Max), true); + VisualizationSceneScalarData::ToogleShading(); } else { diff --git a/lib/vssolution.hpp b/lib/vssolution.hpp index b1cd5e33..29f038c8 100644 --- a/lib/vssolution.hpp +++ b/lib/vssolution.hpp @@ -83,8 +83,6 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData double GetElementLengthScale(int k); public: - Shading shading; - int attr_to_show, bdr_attr_to_show; Array el_attr_to_show, bdr_el_attr_to_show; @@ -166,7 +164,7 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData } virtual void SetShading(Shading, bool); - void ToggleShading(); + virtual void ToggleShading(); void ToggleDrawCP() { draw_cp = !draw_cp; PrepareCP(); } @@ -178,7 +176,6 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData virtual void ToggleAttributes(Array &attr_list); virtual void SetDrawMesh(int i) { drawmesh = i % 3; } - virtual Shading GetShading() { return shading; } virtual int GetDrawMesh() { return drawmesh; } }; diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 508eddf8..ea9ce1d2 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -887,7 +887,7 @@ void VisualizationSceneSolution3d::ToggleShading() { if (GridF) { - SetShading((Shading)(((int)shading+1) % (int)Shading::Max), true); + VisualizationSceneScalarData::ToogleShading(); } else { diff --git a/lib/vssolution3d.hpp b/lib/vssolution3d.hpp index 2179baf2..6f598249 100644 --- a/lib/vssolution3d.hpp +++ b/lib/vssolution3d.hpp @@ -148,9 +148,8 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData // 3 - no arrows (black), 4 - with arrows (black) void ToggleDrawOrdering() { draworder = (draworder+1)%5; } - void ToggleShading(); - Shading GetShading() { return shading; }; virtual void SetShading(Shading, bool); + virtual void ToggleShading(); virtual void SetRefineFactors(int, int); virtual void AutoRefine(); virtual void ToggleAttributes(Array &attr_list); diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index 0149b6c1..7120438f 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -188,7 +188,7 @@ void KeyuPressed() case 0: case 1: if (update && - vsvector->shading == VisualizationSceneSolution::Shading::Noncomforming) + vsvector->GetShading() == VisualizationSceneSolution::Shading::Noncomforming) { vsvector->PrepareVectorField(); SendExposeEvent(); From f4f16cb0da125f37bd789544e9a601e7a2290cb5 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 11:44:44 -0700 Subject: [PATCH 47/57] Added quad mesh lines visualization to 2D scalars. --- glvis.cpp | 5 +- lib/stream_reader.cpp | 2 +- lib/vsdata.cpp | 3 +- lib/vsdata.hpp | 6 +- lib/vssolution.cpp | 146 +++++++++++++++++++++++++++++++++++------- lib/vssolution.hpp | 5 +- lib/vsvector.cpp | 2 +- 7 files changed, 137 insertions(+), 32 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index e48a535f..4091c8d7 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -161,11 +161,12 @@ bool GLVisInitVis(StreamState::FieldType field_type, if (stream_state.normals.Size() > 0) { vs = vss = new VisualizationSceneSolution(*stream_state.mesh, stream_state.sol, - &stream_state.normals); + stream_state.mesh_quad.get(), &stream_state.normals); } else { - vs = vss = new VisualizationSceneSolution(*stream_state.mesh, stream_state.sol); + vs = vss = new VisualizationSceneSolution(*stream_state.mesh, stream_state.sol, + stream_state.mesh_quad.get()); } if (stream_state.grid_f) { diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 4c8e657c..b9022831 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -719,7 +719,7 @@ void StreamState::ResetMeshAndSolution(VisualizationScene* vs) VisualizationSceneSolution *vss = dynamic_cast(vs); grid_f->GetNodalValues(sol); - vss->NewMeshAndSolution(mesh.get(), &sol, grid_f.get()); + vss->NewMeshAndSolution(mesh.get(), mesh_quad.get(), &sol, grid_f.get()); } else { diff --git a/lib/vsdata.cpp b/lib/vsdata.cpp index 1934414f..3c4add44 100644 --- a/lib/vsdata.cpp +++ b/lib/vsdata.cpp @@ -1253,10 +1253,11 @@ void VisualizationSceneScalarData::SetAutoscale(int _autoscale) } VisualizationSceneScalarData::VisualizationSceneScalarData( - Mesh & m, Vector & s) + Mesh & m, Vector & s, Mesh *mc) : a_label_x("x"), a_label_y("y"), a_label_z("z") { mesh = &m; + mesh_coarse = mc; sol = &s; Init(); diff --git a/lib/vsdata.hpp b/lib/vsdata.hpp index d802f48f..92683110 100644 --- a/lib/vsdata.hpp +++ b/lib/vsdata.hpp @@ -68,8 +68,8 @@ class VisualizationSceneScalarData : public VisualizationScene }; protected: - Mesh *mesh; - Vector *sol; + Mesh *mesh{}, *mesh_coarse{}; + Vector *sol{}; double minv, maxv; @@ -151,7 +151,7 @@ class VisualizationSceneScalarData : public VisualizationScene VisualizationSceneScalarData() : a_label_x("x"), a_label_y("y"), a_label_z("z") {} - VisualizationSceneScalarData (Mesh & m, Vector & s); + VisualizationSceneScalarData (Mesh & m, Vector & s, Mesh *mc = NULL); virtual ~VisualizationSceneScalarData(); diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index a3a99270..bab6d0d5 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -423,9 +423,10 @@ VisualizationSceneSolution::VisualizationSceneSolution() } VisualizationSceneSolution::VisualizationSceneSolution( - Mesh &m, Vector &s, Vector *normals) + Mesh &m, Vector &s, Mesh *mc, Vector *normals) { mesh = &m; + mesh_coarse = mc; sol = &s; v_normals = normals; @@ -584,7 +585,7 @@ void VisualizationSceneSolution::ToggleDrawElems() } void VisualizationSceneSolution::NewMeshAndSolution( - Mesh *new_m, Vector *new_sol, GridFunction *new_u) + Mesh *new_m, Mesh *new_mc, Vector *new_sol, GridFunction *new_u) { // If the number of elements changes, recompute the refinement factor if (mesh->GetNE() != new_m->GetNE()) @@ -599,6 +600,7 @@ void VisualizationSceneSolution::NewMeshAndSolution( } } mesh = new_m; + mesh_coarse = new_mc; sol = new_sol; rsol = new_u; @@ -1707,24 +1709,87 @@ void VisualizationSceneSolution::PrepareLines() line_buf.clear(); gl3::GlBuilder lb = line_buf.createBuilder(); - for (i = 0; i < ne; i++) + if (mesh_coarse) { - if (!el_attr_to_show[mesh->GetAttribute(i)-1]) { continue; } + auto &ref = mesh->GetRefinementTransforms(); + IsoparametricTransformation trans; + BiLinear2DFiniteElement fe; + trans.SetFE(&fe); + DenseMatrix emb_pointmat; - lb.glBegin(GL_LINE_LOOP); - mesh->GetPointMatrix (i, pointmat); - mesh->GetElementVertices (i, vertices); - for (j = 0; j < pointmat.Size(); j++) + for (i = 0; i < ne; i++) { - // 1D meshes get rendered flat - double z = GetMinV(); - if (mesh->Dimension() > 1) // In 1D we just put the mesh below the solution + if (!el_attr_to_show[mesh->GetAttribute(i)-1]) { continue; } + + lb.glBegin(GL_LINES); + mesh->GetPointMatrix (i, pointmat); + mesh->GetElementVertices (i, vertices); + + MFEM_ASSERT(pointmat.Size() == 4, "Not a quadrilateral!"); + + //we assume that mesh_course is used only for tensor finite elements, + //like for representation of quadratures, so in 2D it is square + const int geom = + Geometry::Type::SQUARE; //ref.embeddings[i].geom; //<---- bugged!? + const int mat = ref.embeddings[i].matrix; + const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); + trans.SetPointMat(emb_mat); + trans.Transform(fe.GetNodes(), emb_pointmat); + + for (j = 0; j < 4; j++) { - z = LogVal((*sol)(vertices[j])); + int jp1 = (j+1) % 4; + Vector emb_ip1, emb_ip2; + emb_pointmat.GetColumnReference(j, emb_ip1); + emb_pointmat.GetColumnReference(jp1, emb_ip2); + + //check if we are on the parent edge + if (!(( emb_ip1(0) == 0. && emb_ip2(0) == 0.) + || (emb_ip1(0) == 1. && emb_ip2(0) == 1.) + || (emb_ip1(1) == 0. && emb_ip2(1) == 0.) + || (emb_ip1(1) == 1. && emb_ip2(1) == 1.))) + { continue; } + + // 1D meshes get rendered flat + double z1 = GetMinV(); + double z2 = z1; + if (mesh->Dimension() > 1) // In 1D we just put the mesh below the solution + { + z1 = LogVal((*sol)(vertices[j])); + z2 = LogVal((*sol)(vertices[jp1])); + } + + lb.glVertex3d (pointmat(0, j), + pointmat(1, j), + z1); + lb.glVertex3d (pointmat(0, jp1), + pointmat(1, jp1), + z2); } - lb.glVertex3d(pointmat(0, j), pointmat(1, j), z); + lb.glEnd(); + } + } + else + { + for (i = 0; i < ne; i++) + { + if (!el_attr_to_show[mesh->GetAttribute(i)-1]) { continue; } + + lb.glBegin(GL_LINE_LOOP); + mesh->GetPointMatrix (i, pointmat); + mesh->GetElementVertices (i, vertices); + for (j = 0; j < pointmat.Size(); j++) + { + // 1D meshes get rendered flat + double z = GetMinV(); + if (mesh->Dimension() > 1) // In 1D we just put the mesh below the solution + { + z = LogVal((*sol)(vertices[j])); + } + lb.glVertex3d(pointmat(0, j), pointmat(1, j), z); + } + lb.glEnd(); } - lb.glEnd(); } updated_bufs.emplace_back(&line_buf); @@ -2138,14 +2203,51 @@ void VisualizationSceneSolution::PrepareLines3() Array &RE = RefG->RefEdges; lb.glBegin (GL_LINES); - for (k = 0; k < RE.Size()/2; k++) - { - lb.glVertex3d (pointmat(0, RE[2*k]), - pointmat(1, RE[2*k]), - values(RE[2*k])); - lb.glVertex3d (pointmat(0, RE[2*k+1]), - pointmat(1, RE[2*k+1]), - values(RE[2*k+1])); + if (mesh_coarse) + { + auto &ref = mesh->GetRefinementTransforms(); + //we assume that mesh_course is used only for tensor finite elements, + //like for representation of quadratures, so in 2D it is square + const int geom = + Geometry::Type::SQUARE; //ref.embeddings[i].geom; //<---- bugged!? + const int mat = ref.embeddings[i].matrix; + const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); + IsoparametricTransformation trans; + BiLinear2DFiniteElement fe; + trans.SetFE(&fe); + trans.SetPointMat(emb_mat); + for (k = 0; k < RE.Size()/2; k++) + { + Vector emb_ip1, emb_ip2; + trans.Transform(RefG->RefPts[RE[2*k]], emb_ip1); + trans.Transform(RefG->RefPts[RE[2*k+1]], emb_ip2); + + //check if we are on the parent edge + if (!(( emb_ip1(0) == 0. && emb_ip2(0) == 0.) + || (emb_ip1(0) == 1. && emb_ip2(0) == 1.) + || (emb_ip1(1) == 0. && emb_ip2(1) == 0.) + || (emb_ip1(1) == 1. && emb_ip2(1) == 1.))) + { continue; } + + lb.glVertex3d (pointmat(0, RE[2*k]), + pointmat(1, RE[2*k]), + values(RE[2*k])); + lb.glVertex3d (pointmat(0, RE[2*k+1]), + pointmat(1, RE[2*k+1]), + values(RE[2*k+1])); + } + } + else + { + for (k = 0; k < RE.Size()/2; k++) + { + lb.glVertex3d (pointmat(0, RE[2*k]), + pointmat(1, RE[2*k]), + values(RE[2*k])); + lb.glVertex3d (pointmat(0, RE[2*k+1]), + pointmat(1, RE[2*k+1]), + values(RE[2*k+1])); + } } lb.glEnd(); } diff --git a/lib/vssolution.hpp b/lib/vssolution.hpp index 29f038c8..d46c6ef0 100644 --- a/lib/vssolution.hpp +++ b/lib/vssolution.hpp @@ -87,7 +87,8 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData Array el_attr_to_show, bdr_el_attr_to_show; VisualizationSceneSolution(); - VisualizationSceneSolution(Mesh &m, Vector &s, Vector *normals = NULL); + VisualizationSceneSolution(Mesh &m, Vector &s, Mesh *mc = NULL, + Vector *normals = NULL); virtual ~VisualizationSceneSolution(); @@ -95,7 +96,7 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData void SetGridFunction(GridFunction & u) { rsol = &u; } - void NewMeshAndSolution(Mesh *new_m, Vector *new_sol, + void NewMeshAndSolution(Mesh *new_m, Mesh *new_mc, Vector *new_sol, GridFunction *new_u = NULL); virtual void SetNewScalingFromBox(); diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index 7120438f..66ae540d 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -459,7 +459,7 @@ void VisualizationSceneVector::NewMeshAndSolution(GridFunction &vgf) (*sol)(i) = Vec2Scalar((*solx)(i), (*soly)(i)); } - VisualizationSceneSolution::NewMeshAndSolution(mesh, sol, &vgf); + VisualizationSceneSolution::NewMeshAndSolution(mesh, NULL, sol, &vgf); if (autoscale) { From 4cb452c3fd092d98b3a0a940e9650b1a438d5f19 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 11:59:35 -0700 Subject: [PATCH 48/57] Added quad mesh lines visualization to 2D vectors. --- glvis.cpp | 2 +- lib/stream_reader.cpp | 2 +- lib/vsvector.cpp | 8 +++++--- lib/vsvector.hpp | 4 ++-- 4 files changed, 9 insertions(+), 7 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index 4091c8d7..48e92284 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -232,7 +232,7 @@ bool GLVisInitVis(StreamState::FieldType field_type, else { vs = new VisualizationSceneVector(*stream_state.mesh, stream_state.solu, - stream_state.solv); + stream_state.solv, stream_state.mesh_quad.get()); } } else if (stream_state.mesh->SpaceDimension() == 3) diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index b9022831..67fdf5b9 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -725,7 +725,7 @@ void StreamState::ResetMeshAndSolution(VisualizationScene* vs) { VisualizationSceneVector *vsv = dynamic_cast(vs); - vsv->NewMeshAndSolution(*grid_f); + vsv->NewMeshAndSolution(*grid_f, mesh_quad.get()); } } else diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index 66ae540d..96e3b317 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -255,9 +255,10 @@ const char *Vec2ScalarNames[7] = }; VisualizationSceneVector::VisualizationSceneVector(Mesh & m, - Vector & sx, Vector & sy) + Vector & sx, Vector & sy, Mesh *mc) { mesh = &m; + mesh_coarse = mc; solx = &sx; soly = &sy; @@ -405,7 +406,7 @@ void VisualizationSceneVector::CycleVec2Scalar(int print) } } -void VisualizationSceneVector::NewMeshAndSolution(GridFunction &vgf) +void VisualizationSceneVector::NewMeshAndSolution(GridFunction &vgf, Mesh *mc) { delete sol; @@ -436,6 +437,7 @@ void VisualizationSceneVector::NewMeshAndSolution(GridFunction &vgf) } } mesh = new_mesh; + mesh_coarse = mc; solx = new Vector(mesh->GetNV()); soly = new Vector(mesh->GetNV()); @@ -459,7 +461,7 @@ void VisualizationSceneVector::NewMeshAndSolution(GridFunction &vgf) (*sol)(i) = Vec2Scalar((*solx)(i), (*soly)(i)); } - VisualizationSceneSolution::NewMeshAndSolution(mesh, NULL, sol, &vgf); + VisualizationSceneSolution::NewMeshAndSolution(mesh, mesh_coarse, sol, &vgf); if (autoscale) { diff --git a/lib/vsvector.hpp b/lib/vsvector.hpp index ea41222a..5593f466 100644 --- a/lib/vsvector.hpp +++ b/lib/vsvector.hpp @@ -46,10 +46,10 @@ class VisualizationSceneVector : public VisualizationSceneSolution IsoparametricTransformation T0; public: - VisualizationSceneVector(Mesh &m, Vector &sx, Vector &sy); + VisualizationSceneVector(Mesh &m, Vector &sx, Vector &sy, Mesh *mc = NULL); VisualizationSceneVector(GridFunction &vgf); - void NewMeshAndSolution(GridFunction &vgf); + void NewMeshAndSolution(GridFunction &vgf, Mesh *mc = NULL); virtual ~VisualizationSceneVector(); From 302802dc178e4a877e8e46f1fc0a019e64830ecb Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 14:40:27 -0700 Subject: [PATCH 49/57] Fixed move of shading. --- lib/vssolution3d.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/lib/vssolution3d.hpp b/lib/vssolution3d.hpp index 6f598249..b1889516 100644 --- a/lib/vssolution3d.hpp +++ b/lib/vssolution3d.hpp @@ -23,7 +23,6 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData protected: int drawmesh, drawelems, draworder; - Shading shading; int cplane; int cp_drawmesh, cp_drawelems, drawlsurf; // Algorithm used to draw the cutting plane when shading is 2 and cplane is 1 From 055433735ccade2f5e8679a2a15a7f5e76fb359d Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 17:24:11 -0700 Subject: [PATCH 50/57] Added quad lines for scalar 3D. --- glvis.cpp | 2 +- lib/stream_reader.cpp | 2 +- lib/vssolution3d.cpp | 159 +++++++++++++++++++++++++++++++++++++++--- lib/vssolution3d.hpp | 4 +- 4 files changed, 154 insertions(+), 13 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index 48e92284..e87fc101 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -186,7 +186,7 @@ bool GLVisInitVis(StreamState::FieldType field_type, { VisualizationSceneSolution3d * vss; vs = vss = new VisualizationSceneSolution3d(*stream_state.mesh, - stream_state.sol); + stream_state.sol, stream_state.mesh_quad.get()); if (stream_state.grid_f) { vss->SetGridFunction(stream_state.grid_f.get()); diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 67fdf5b9..8bc8c475 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -735,7 +735,7 @@ void StreamState::ResetMeshAndSolution(VisualizationScene* vs) VisualizationSceneSolution3d *vss = dynamic_cast(vs); grid_f->GetNodalValues(sol); - vss->NewMeshAndSolution(mesh.get(), &sol, grid_f.get()); + vss->NewMeshAndSolution(mesh.get(), mesh_quad.get(), &sol, grid_f.get()); } else { diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index ea9ce1d2..186e10d9 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -688,9 +688,11 @@ static void KeyF10Pressed() VisualizationSceneSolution3d::VisualizationSceneSolution3d() {} -VisualizationSceneSolution3d::VisualizationSceneSolution3d(Mesh &m, Vector &s) +VisualizationSceneSolution3d::VisualizationSceneSolution3d(Mesh &m, Vector &s, + Mesh *mc) { mesh = &m; + mesh_coarse = mc; sol = &s; GridF = NULL; @@ -818,7 +820,7 @@ VisualizationSceneSolution3d::~VisualizationSceneSolution3d() } void VisualizationSceneSolution3d::NewMeshAndSolution( - Mesh *new_m, Vector *new_sol, GridFunction *new_u) + Mesh *new_m, Mesh *new_mc, Vector *new_sol, GridFunction *new_u) { if (mesh->GetNV() != new_m->GetNV()) { @@ -839,6 +841,7 @@ void VisualizationSceneSolution3d::NewMeshAndSolution( } } mesh = new_m; + mesh_coarse = new_mc; sol = new_sol; GridF = new_u; FindNodePos(); @@ -2070,15 +2073,88 @@ void VisualizationSceneSolution3d::PrepareLines() switch (drawmesh) { case 1: - line.glBegin(GL_LINE_LOOP); + { + if (mesh_coarse) + { + int f, o, e1; + FaceElementTransformations *ftr; + if (dim == 3) + { + mesh->GetBdrElementFace(i, &f, &o); + ftr = mesh->GetFaceElementTransformations(f); + e1 = ftr->Elem1No; + } + else + { + ftr = NULL; + e1 = i; + } + auto &ref = mesh->GetRefinementTransforms(); + IsoparametricTransformation trans; + BiLinear2DFiniteElement fe_face; + TriLinear3DFiniteElement fe; + trans.SetFE(&fe); + DenseMatrix emb_pointmat; + + MFEM_ASSERT(pointmat.Size() == 4, "Not a quadrilateral!"); + + //we assume that mesh_course is used only for tensor finite elements, + //like for representation of quadratures, so in 2D it is square + const int geom = (dim == 3)?(Geometry::Type::CUBE) + :(Geometry::Type::SQUARE); //ref.embeddings[e1].geom; //<---- bugged!? + const int mat = ref.embeddings[e1].matrix; + const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); + trans.SetPointMat(emb_mat); + if (dim == 3) + { + IntegrationRule nodes3d(4); + ftr->Loc1.Transform(fe_face.GetNodes(), nodes3d); + trans.Transform(nodes3d, emb_pointmat); + } + else + { + trans.Transform(fe_face.GetNodes(), emb_pointmat); + } - for (j = 0; j < pointmat.Size(); j++) + line.glBegin(GL_LINES); + + for (j = 0; j < 4; j++) + { + int jp1 = (j+1) % 4; + Vector emb_ip1, emb_ip2; + emb_pointmat.GetColumnReference(j, emb_ip1); + emb_pointmat.GetColumnReference(jp1, emb_ip2); + + //check if we are on the outer edge + int inter = 0; + for (int d = 0; d < dim; d++) + if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) + || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) + { inter++; } + if (inter != 1) { continue; } + + line.glVertex3d (pointmat(0, j), + pointmat(1, j), + pointmat(2, j)); + line.glVertex3d (pointmat(0, jp1), + pointmat(1, jp1), + pointmat(2, jp1)); + } + + line.glEnd(); + } + else { - line.glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j)); + line.glBegin(GL_LINE_LOOP); + + for (j = 0; j < pointmat.Size(); j++) + { + line.glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j)); + } + line.glEnd(); } - line.glEnd(); break; - + } case 2: for (j = 0; j < pointmat.Size(); j++) { @@ -2220,9 +2296,74 @@ void VisualizationSceneSolution3d::PrepareLines2() Array &REdges = RefG->RefEdges; line.glBegin(GL_LINES); - for (int k = 0; k < REdges.Size(); k++) + + if (mesh_coarse) + { + int f, o, e1; + FaceElementTransformations *ftr; + if (dim == 3) + { + mesh->GetBdrElementFace(i, &f, &o); + ftr = mesh->GetFaceElementTransformations(f); + e1 = ftr->Elem1No; + } + else + { + ftr = NULL; + e1 = i; + } + auto &ref = mesh->GetRefinementTransforms(); + IsoparametricTransformation trans; + TriLinear3DFiniteElement fe; + trans.SetFE(&fe); + + //we assume that mesh_course is used only for tensor finite elements, + //like for representation of quadratures, so in 2D it is square + const int geom = (dim == 3)?(Geometry::Type::CUBE) + :(Geometry::Type::SQUARE); //ref.embeddings[e1].geom; //<---- bugged!? + const int mat = ref.embeddings[e1].matrix; + const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); + trans.SetPointMat(emb_mat); + + line.glBegin(GL_LINES); + + for (int k = 0; k < REdges.Size()/2; k++) + { + Vector emb_ip1, emb_ip2; + if (dim == 3) + { + IntegrationPoint ip1_3d, ip2_3d; + ftr->Loc1.Transform(RefG->RefPts[REdges[2*k]], ip1_3d); + ftr->Loc1.Transform(RefG->RefPts[REdges[2*k+1]], ip2_3d); + trans.Transform(ip1_3d, emb_ip1); + trans.Transform(ip2_3d, emb_ip2); + } + else + { + trans.Transform(RefG->RefPts[REdges[2*k]], emb_ip1); + trans.Transform(RefG->RefPts[REdges[2*k+1]], emb_ip2); + } + + //check if we are on the outer edge + int inter = 0; + for (int d = 0; d < dim; d++) + if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) + || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) + { inter++; } + if (inter != 1) { continue; } + + line.glVertex3dv(&pointmat(0, REdges[2*k])); + line.glVertex3dv(&pointmat(0, REdges[2*k+1])); + } + + line.glEnd(); + } + else { - line.glVertex3dv(&pointmat(0, REdges[k])); + for (int k = 0; k < REdges.Size(); k++) + { + line.glVertex3dv(&pointmat(0, REdges[k])); + } } line.glEnd(); } diff --git a/lib/vssolution3d.hpp b/lib/vssolution3d.hpp index b1889516..205b109e 100644 --- a/lib/vssolution3d.hpp +++ b/lib/vssolution3d.hpp @@ -112,11 +112,11 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData Array bdr_attr_to_show; VisualizationSceneSolution3d(); - VisualizationSceneSolution3d(Mesh & m, Vector & s); + VisualizationSceneSolution3d(Mesh & m, Vector & s, Mesh *mc); void SetGridFunction (GridFunction *gf) { GridF = gf; } - void NewMeshAndSolution(Mesh *new_m, Vector *new_sol, + void NewMeshAndSolution(Mesh *new_m, Mesh *new_mc, Vector *new_sol, GridFunction *new_u = NULL); virtual ~VisualizationSceneSolution3d(); From 02e8f73a15aafd5e69eb2e4b9deda732c529ad26 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 17:52:24 -0700 Subject: [PATCH 51/57] Added quad lines for scalar 3D cplane=2. --- lib/vssolution3d.cpp | 116 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 110 insertions(+), 6 deletions(-) diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 186e10d9..55b5e099 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -3246,14 +3246,64 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() switch (cp_drawmesh) { case 1: - // glBegin(GL_POLYGON); - line.glBegin(GL_LINE_LOOP); - for (j = 0; j < nodes.Size(); j++) + { + if (mesh_coarse) + { + FaceElementTransformations *ftr; + ftr = mesh->GetFaceElementTransformations(i); + auto &ref = mesh->GetRefinementTransforms(); + IsoparametricTransformation trans; + BiLinear2DFiniteElement fe_face; + TriLinear3DFiniteElement fe; + trans.SetFE(&fe); + DenseMatrix emb_pointmat; + + //we assume that mesh_course is used only for tensor finite elements, + //like for representation of quadratures, so in 2D it is square + const int geom = + Geometry::Type::CUBE; //ref.embeddings[e1].geom; //<---- bugged!? + const int mat = ref.embeddings[e1].matrix; + const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); + trans.SetPointMat(emb_mat); + IntegrationRule nodes3d(4); + ftr->Loc1.Transform(fe_face.GetNodes(), nodes3d); + trans.Transform(nodes3d, emb_pointmat); + + line.glBegin(GL_LINES); + + for (j = 0; j < 4; j++) + { + int jp1 = (j+1) % 4; + Vector emb_ip1, emb_ip2; + emb_pointmat.GetColumnReference(j, emb_ip1); + emb_pointmat.GetColumnReference(jp1, emb_ip2); + + //check if we are on the outer edge + int inter = 0; + for (int d = 0; d < 3; d++) + if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) + || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) + { inter++; } + if (inter != 1) { continue; } + + line.glVertex3dv(point[j]); + line.glVertex3dv(point[jp1]); + } + + line.glEnd(); + } + else { - line.glVertex3dv (point[j]); + // glBegin(GL_POLYGON); + line.glBegin(GL_LINE_LOOP); + for (j = 0; j < nodes.Size(); j++) + { + line.glVertex3dv (point[j]); + } + line.glEnd(); } - line.glEnd(); break; + } case 2: DrawPolygonLevelLines(line, point[0], nodes.Size(), level, false); break; @@ -3278,8 +3328,62 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() switch (cp_drawmesh) { case 1: - DrawRefinedSurfEdges (n, pointmat, values, RefG->RefEdges); + { + if (mesh_coarse) + { + FaceElementTransformations *ftr; + ftr = mesh->GetFaceElementTransformations(i); + auto &ref = mesh->GetRefinementTransforms(); + IsoparametricTransformation trans; + BiLinear2DFiniteElement fe_face; + TriLinear3DFiniteElement fe; + trans.SetFE(&fe); + DenseMatrix emb_pointmat; + + //we assume that mesh_course is used only for tensor finite elements, + //like for representation of quadratures, so in 2D it is square + const int geom = + Geometry::Type::CUBE; //ref.embeddings[e1].geom; //<---- bugged!? + const int mat = ref.embeddings[e1].matrix; + const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); + trans.SetPointMat(emb_mat); + IntegrationRule nodes3d(4); + ftr->Loc1.Transform(fe_face.GetNodes(), nodes3d); + trans.Transform(nodes3d, emb_pointmat); + + gl3::GlBuilder line = cplines_buf.createBuilder(); + line.glBegin(GL_LINES); + + auto &RE = RefG->RefEdges; + for (int k = 0; k < RE.Size()/2; k++) + { + IntegrationPoint ip1_3d, ip2_3d; + Vector emb_ip1, emb_ip2; + ftr->Loc1.Transform(RefG->RefPts[RE[2*k]], ip1_3d); + ftr->Loc1.Transform(RefG->RefPts[RE[2*k+1]], ip2_3d); + trans.Transform(ip1_3d, emb_ip1); + trans.Transform(ip2_3d, emb_ip2); + + //check if we are on the outer edge + int inter = 0; + for (int d = 0; d < 3; d++) + if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) + || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) + { inter++; } + if (inter != 1) { continue; } + + line.glVertex3dv(&pointmat(0, RE[2*k])); + line.glVertex3dv(&pointmat(0, RE[2*k+1])); + } + + line.glEnd(); + } + else + { + DrawRefinedSurfEdges (n, pointmat, values, RefG->RefEdges); + } break; + } case 2: DrawRefinedSurfLevelLines (n, pointmat, values, RefG->RefGeoms); From 04938f067c87a2f40ad0cf544b986c465118e99d Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 19:11:24 -0700 Subject: [PATCH 52/57] Fixed cutting plane outer edge detection. --- lib/vssolution3d.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 55b5e099..58334d20 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -3284,6 +3284,8 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) { inter++; } + if (ref.embeddings[e1].parent == ref.embeddings[e2].parent) + { inter--; } if (inter != 1) { continue; } line.glVertex3dv(point[j]); @@ -3370,6 +3372,8 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) { inter++; } + if (ref.embeddings[e1].parent == ref.embeddings[e2].parent) + { inter--; } if (inter != 1) { continue; } line.glVertex3dv(&pointmat(0, RE[2*k])); From 55694c85c66ca4c79234399a884eab4255b6cc96 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 22:11:48 -0700 Subject: [PATCH 53/57] Unified codebase in 3D scalars. --- lib/vssolution3d.cpp | 365 ++++++++++++++++--------------------------- lib/vssolution3d.hpp | 8 + 2 files changed, 146 insertions(+), 227 deletions(-) diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 58334d20..6d275798 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -1410,6 +1410,128 @@ void VisualizationSceneSolution3d::DrawRefinedSurfEdges( line.glEnd(); } +void VisualizationSceneSolution3d::DrawBdrElCoarseSurfEdges( + gl3::GlBuilder &line, int be, DenseMatrix &pointmat, const IntegrationRule *ir, + Array *idxs) +{ + const int dim = mesh->Dimension(); + int f, e1; + if (dim == 3) + { + int o; + mesh->GetBdrElementFace(be, &f, &o); + e1 = -1; + } + else + { + f = -1; + e1 = be; + } + DrawCoarseSurfEdges(line, f, e1, -1, pointmat, ir, idxs); +} + +void VisualizationSceneSolution3d::DrawFaceCoarseSurfEdges( + gl3::GlBuilder &line, int f, DenseMatrix &pointmat, const IntegrationRule *ir, + Array *idxs) +{ + DrawCoarseSurfEdges(line, f, -1, -1, pointmat, ir, idxs); +} + +void VisualizationSceneSolution3d::DrawCoarseSurfEdges( + gl3::GlBuilder &line, int f, int e1, int e2, DenseMatrix &pointmat, + const IntegrationRule *ir, + Array *idxs) +{ + MFEM_ASSERT(mesh_coarse, "Cannot be used without the coarse mesh"); + MFEM_ASSERT(mesh->GetLastOperation() == Mesh::Operation::REFINE, + "Not a refined mesh"); + MFEM_ASSERT(!idxs || idxs->Size() == 8, "Not expected 4 edges"); + + const int dim = mesh->Dimension(); + FaceElementTransformations *ftr; + if (f >= 0) + { + ftr = mesh->GetFaceElementTransformations(f); + e1 = ftr->Elem1No; + e2 = ftr->Elem2No; + } + auto &ref = mesh->GetRefinementTransforms(); + IsoparametricTransformation trans; + BiLinear2DFiniteElement fe_face; + TriLinear3DFiniteElement fe; + trans.SetFE(&fe); + DenseMatrix emb_pointmat; + + //we assume that mesh_course is used only for tensor finite elements, + //like for representation of quadratures, so in 2D it is square + const int geom = (dim == 3)?(Geometry::Type::CUBE) + :(Geometry::Type::SQUARE); //ref.embeddings[e1].geom; //<---- bugged!? + const int mat = ref.embeddings[e1].matrix; + const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); + trans.SetPointMat(emb_mat); + if (!ir) + { + if (dim == 3) + { + IntegrationRule nodes3d(4); + ftr->Loc1.Transform(fe_face.GetNodes(), nodes3d); + trans.Transform(nodes3d, emb_pointmat); + } + else + { + trans.Transform(fe_face.GetNodes(), emb_pointmat); + } + } + + line.glBegin(GL_LINES); + + for (int k = 0; k < 4; k++) + { + int j, jp1; + Vector emb_ip1, emb_ip2; + if (ir && idxs) + { + j = (*idxs)[2*k]; + jp1 = (*idxs)[2*k+1]; + if (dim == 3) + { + IntegrationPoint ip1_3d, ip2_3d; + ftr->Loc1.Transform((*ir)[j], ip1_3d); + ftr->Loc1.Transform((*ir)[jp1], ip2_3d); + trans.Transform(ip1_3d, emb_ip1); + trans.Transform(ip2_3d, emb_ip2); + } + else + { + trans.Transform((*ir)[j], emb_ip1); + trans.Transform((*ir)[jp1], emb_ip2); + } + } + else + { + j = k; + jp1 = (k+1) % 4; + emb_pointmat.GetColumnReference(j, emb_ip1); + emb_pointmat.GetColumnReference(jp1, emb_ip2); + } + + //check if we are on the outer edge + int inter = 0; + for (int d = 0; d < 3; d++) + if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) + || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) + { inter++; } + if (e2>= 0 && ref.embeddings[e1].parent == ref.embeddings[e2].parent) + { inter--; } + if (inter != 1) { continue; } + + line.glVertex3dv(&pointmat(0, j)); + line.glVertex3dv(&pointmat(0, jp1)); + } + + line.glEnd(); +} + void VisualizationSceneSolution3d::DrawRefinedSurfLevelLines( int n, DenseMatrix &pointmat, Vector &values, Array &RefGeoms) { @@ -2024,14 +2146,13 @@ void VisualizationSceneSolution3d::PrepareLines() int dim = mesh->Dimension(); int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE(); - int i, j, k; DenseMatrix pointmat; line_buf.clear(); Array vertices; - for (i = 0; i < ne; i++) + for (int i = 0; i < ne; i++) { if (dim == 3) { @@ -2076,78 +2197,13 @@ void VisualizationSceneSolution3d::PrepareLines() { if (mesh_coarse) { - int f, o, e1; - FaceElementTransformations *ftr; - if (dim == 3) - { - mesh->GetBdrElementFace(i, &f, &o); - ftr = mesh->GetFaceElementTransformations(f); - e1 = ftr->Elem1No; - } - else - { - ftr = NULL; - e1 = i; - } - auto &ref = mesh->GetRefinementTransforms(); - IsoparametricTransformation trans; - BiLinear2DFiniteElement fe_face; - TriLinear3DFiniteElement fe; - trans.SetFE(&fe); - DenseMatrix emb_pointmat; - - MFEM_ASSERT(pointmat.Size() == 4, "Not a quadrilateral!"); - - //we assume that mesh_course is used only for tensor finite elements, - //like for representation of quadratures, so in 2D it is square - const int geom = (dim == 3)?(Geometry::Type::CUBE) - :(Geometry::Type::SQUARE); //ref.embeddings[e1].geom; //<---- bugged!? - const int mat = ref.embeddings[e1].matrix; - const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); - trans.SetPointMat(emb_mat); - if (dim == 3) - { - IntegrationRule nodes3d(4); - ftr->Loc1.Transform(fe_face.GetNodes(), nodes3d); - trans.Transform(nodes3d, emb_pointmat); - } - else - { - trans.Transform(fe_face.GetNodes(), emb_pointmat); - } - - line.glBegin(GL_LINES); - - for (j = 0; j < 4; j++) - { - int jp1 = (j+1) % 4; - Vector emb_ip1, emb_ip2; - emb_pointmat.GetColumnReference(j, emb_ip1); - emb_pointmat.GetColumnReference(jp1, emb_ip2); - - //check if we are on the outer edge - int inter = 0; - for (int d = 0; d < dim; d++) - if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) - || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) - { inter++; } - if (inter != 1) { continue; } - - line.glVertex3d (pointmat(0, j), - pointmat(1, j), - pointmat(2, j)); - line.glVertex3d (pointmat(0, jp1), - pointmat(1, jp1), - pointmat(2, jp1)); - } - - line.glEnd(); + DrawBdrElCoarseSurfEdges(line, i, pointmat); } else { line.glBegin(GL_LINE_LOOP); - for (j = 0; j < pointmat.Size(); j++) + for (int j = 0; j < pointmat.Size(); j++) { line.glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j)); } @@ -2156,9 +2212,9 @@ void VisualizationSceneSolution3d::PrepareLines() break; } case 2: - for (j = 0; j < pointmat.Size(); j++) + for (int j = 0; j < pointmat.Size(); j++) { - for (k = 0; k < 3; k++) + for (int k = 0; k < 3; k++) { point[j][k] = pointmat(k,j); } @@ -2299,64 +2355,7 @@ void VisualizationSceneSolution3d::PrepareLines2() if (mesh_coarse) { - int f, o, e1; - FaceElementTransformations *ftr; - if (dim == 3) - { - mesh->GetBdrElementFace(i, &f, &o); - ftr = mesh->GetFaceElementTransformations(f); - e1 = ftr->Elem1No; - } - else - { - ftr = NULL; - e1 = i; - } - auto &ref = mesh->GetRefinementTransforms(); - IsoparametricTransformation trans; - TriLinear3DFiniteElement fe; - trans.SetFE(&fe); - - //we assume that mesh_course is used only for tensor finite elements, - //like for representation of quadratures, so in 2D it is square - const int geom = (dim == 3)?(Geometry::Type::CUBE) - :(Geometry::Type::SQUARE); //ref.embeddings[e1].geom; //<---- bugged!? - const int mat = ref.embeddings[e1].matrix; - const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); - trans.SetPointMat(emb_mat); - - line.glBegin(GL_LINES); - - for (int k = 0; k < REdges.Size()/2; k++) - { - Vector emb_ip1, emb_ip2; - if (dim == 3) - { - IntegrationPoint ip1_3d, ip2_3d; - ftr->Loc1.Transform(RefG->RefPts[REdges[2*k]], ip1_3d); - ftr->Loc1.Transform(RefG->RefPts[REdges[2*k+1]], ip2_3d); - trans.Transform(ip1_3d, emb_ip1); - trans.Transform(ip2_3d, emb_ip2); - } - else - { - trans.Transform(RefG->RefPts[REdges[2*k]], emb_ip1); - trans.Transform(RefG->RefPts[REdges[2*k+1]], emb_ip2); - } - - //check if we are on the outer edge - int inter = 0; - for (int d = 0; d < dim; d++) - if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) - || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) - { inter++; } - if (inter != 1) { continue; } - - line.glVertex3dv(&pointmat(0, REdges[2*k])); - line.glVertex3dv(&pointmat(0, REdges[2*k+1])); - } - - line.glEnd(); + DrawBdrElCoarseSurfEdges(line, i, pointmat, &RefG->RefPts, &REdges); } else { @@ -3207,7 +3206,7 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines() void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() { int i, j, n = 0; - double point[4][4], *coord; + double *coord; DenseMatrix pointmat; Vector values; RefinedGeometry *RefG; @@ -3234,13 +3233,14 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() if (shading != Shading::Noncomforming) { mesh -> GetFaceVertices (i, nodes); + pointmat.SetSize(4, nodes.Size()); for (j = 0; j < nodes.Size(); j++) { coord = mesh -> GetVertex(nodes[j]); - point[j][0] = coord[0]; - point[j][1] = coord[1]; - point[j][2] = coord[2]; - point[j][3] = (*sol)(nodes[j]); + pointmat(0, j) = coord[0]; + pointmat(1, j) = coord[1]; + pointmat(2, j) = coord[2]; + pointmat(3, j) = (*sol)(nodes[j]); } gl3::GlBuilder line = cplines_buf.createBuilder(); switch (cp_drawmesh) @@ -3249,50 +3249,7 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() { if (mesh_coarse) { - FaceElementTransformations *ftr; - ftr = mesh->GetFaceElementTransformations(i); - auto &ref = mesh->GetRefinementTransforms(); - IsoparametricTransformation trans; - BiLinear2DFiniteElement fe_face; - TriLinear3DFiniteElement fe; - trans.SetFE(&fe); - DenseMatrix emb_pointmat; - - //we assume that mesh_course is used only for tensor finite elements, - //like for representation of quadratures, so in 2D it is square - const int geom = - Geometry::Type::CUBE; //ref.embeddings[e1].geom; //<---- bugged!? - const int mat = ref.embeddings[e1].matrix; - const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); - trans.SetPointMat(emb_mat); - IntegrationRule nodes3d(4); - ftr->Loc1.Transform(fe_face.GetNodes(), nodes3d); - trans.Transform(nodes3d, emb_pointmat); - - line.glBegin(GL_LINES); - - for (j = 0; j < 4; j++) - { - int jp1 = (j+1) % 4; - Vector emb_ip1, emb_ip2; - emb_pointmat.GetColumnReference(j, emb_ip1); - emb_pointmat.GetColumnReference(jp1, emb_ip2); - - //check if we are on the outer edge - int inter = 0; - for (int d = 0; d < 3; d++) - if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) - || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) - { inter++; } - if (ref.embeddings[e1].parent == ref.embeddings[e2].parent) - { inter--; } - if (inter != 1) { continue; } - - line.glVertex3dv(point[j]); - line.glVertex3dv(point[jp1]); - } - - line.glEnd(); + DrawFaceCoarseSurfEdges(line, i, pointmat); } else { @@ -3300,14 +3257,14 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() line.glBegin(GL_LINE_LOOP); for (j = 0; j < nodes.Size(); j++) { - line.glVertex3dv (point[j]); + line.glVertex3dv(&pointmat(0,j)); } line.glEnd(); } break; } case 2: - DrawPolygonLevelLines(line, point[0], nodes.Size(), level, false); + DrawPolygonLevelLines(line, pointmat.GetData(), nodes.Size(), level, false); break; } } @@ -3333,54 +3290,8 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() { if (mesh_coarse) { - FaceElementTransformations *ftr; - ftr = mesh->GetFaceElementTransformations(i); - auto &ref = mesh->GetRefinementTransforms(); - IsoparametricTransformation trans; - BiLinear2DFiniteElement fe_face; - TriLinear3DFiniteElement fe; - trans.SetFE(&fe); - DenseMatrix emb_pointmat; - - //we assume that mesh_course is used only for tensor finite elements, - //like for representation of quadratures, so in 2D it is square - const int geom = - Geometry::Type::CUBE; //ref.embeddings[e1].geom; //<---- bugged!? - const int mat = ref.embeddings[e1].matrix; - const DenseMatrix &emb_mat = ref.point_matrices[geom](mat); - trans.SetPointMat(emb_mat); - IntegrationRule nodes3d(4); - ftr->Loc1.Transform(fe_face.GetNodes(), nodes3d); - trans.Transform(nodes3d, emb_pointmat); - gl3::GlBuilder line = cplines_buf.createBuilder(); - line.glBegin(GL_LINES); - - auto &RE = RefG->RefEdges; - for (int k = 0; k < RE.Size()/2; k++) - { - IntegrationPoint ip1_3d, ip2_3d; - Vector emb_ip1, emb_ip2; - ftr->Loc1.Transform(RefG->RefPts[RE[2*k]], ip1_3d); - ftr->Loc1.Transform(RefG->RefPts[RE[2*k+1]], ip2_3d); - trans.Transform(ip1_3d, emb_ip1); - trans.Transform(ip2_3d, emb_ip2); - - //check if we are on the outer edge - int inter = 0; - for (int d = 0; d < 3; d++) - if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) - || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) - { inter++; } - if (ref.embeddings[e1].parent == ref.embeddings[e2].parent) - { inter--; } - if (inter != 1) { continue; } - - line.glVertex3dv(&pointmat(0, RE[2*k])); - line.glVertex3dv(&pointmat(0, RE[2*k+1])); - } - - line.glEnd(); + DrawFaceCoarseSurfEdges(line, i, pointmat, &RefG->RefPts, &RefG->RefEdges); } else { diff --git a/lib/vssolution3d.hpp b/lib/vssolution3d.hpp index 205b109e..1b3336ce 100644 --- a/lib/vssolution3d.hpp +++ b/lib/vssolution3d.hpp @@ -59,6 +59,14 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData void DrawRefinedSurfEdges (int n, DenseMatrix &pointmat, Vector &values, Array &RefEdges, int part = -1); + void DrawBdrElCoarseSurfEdges(gl3::GlBuilder &line, int be, + DenseMatrix &pointmat, const IntegrationRule *ir = NULL, + Array *idxs = NULL); + void DrawFaceCoarseSurfEdges(gl3::GlBuilder &line, int f, DenseMatrix &pointmat, + const IntegrationRule *ir = NULL, Array *idxs = NULL); + void DrawCoarseSurfEdges(gl3::GlBuilder &line, int f, int e1, int e2, + DenseMatrix &pointmat, const IntegrationRule *ir = NULL, + Array *idxs = NULL); void LiftRefinedSurf (int n, DenseMatrix &pointmat, Vector &values, int *RG); void DrawTetLevelSurf(gl3::GlDrawable& target, const DenseMatrix &verts, From 1157a975b930c9b35648683eed3b3da461edcea8 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 22:47:04 -0700 Subject: [PATCH 54/57] Fixed 3D coarse meshlines for refined edges. --- lib/vssolution3d.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 6d275798..b045a530 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -1445,7 +1445,6 @@ void VisualizationSceneSolution3d::DrawCoarseSurfEdges( MFEM_ASSERT(mesh_coarse, "Cannot be used without the coarse mesh"); MFEM_ASSERT(mesh->GetLastOperation() == Mesh::Operation::REFINE, "Not a refined mesh"); - MFEM_ASSERT(!idxs || idxs->Size() == 8, "Not expected 4 edges"); const int dim = mesh->Dimension(); FaceElementTransformations *ftr; @@ -1485,7 +1484,9 @@ void VisualizationSceneSolution3d::DrawCoarseSurfEdges( line.glBegin(GL_LINES); - for (int k = 0; k < 4; k++) + const int k_max = (idxs)?(idxs->Size()/2):(4); + + for (int k = 0; k < k_max; k++) { int j, jp1; Vector emb_ip1, emb_ip2; From 5d273849bd61b28294fd4270b1db755af2904aba Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 18 Jul 2024 22:51:25 -0700 Subject: [PATCH 55/57] Added quad mesh lines to vector 3D. --- glvis.cpp | 6 ++- lib/stream_reader.cpp | 2 +- lib/vsvector3d.cpp | 90 ++++++++++++++++++++++++++++++------------- lib/vsvector3d.hpp | 7 ++-- 4 files changed, 72 insertions(+), 33 deletions(-) diff --git a/glvis.cpp b/glvis.cpp index e87fc101..82700368 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -240,12 +240,14 @@ bool GLVisInitVis(StreamState::FieldType field_type, if (stream_state.grid_f) { stream_state.ProjectVectorFEGridFunction(); - vs = new VisualizationSceneVector3d(*stream_state.grid_f); + vs = new VisualizationSceneVector3d(*stream_state.grid_f, + stream_state.mesh_quad.get()); } else { vs = new VisualizationSceneVector3d(*stream_state.mesh, stream_state.solu, - stream_state.solv, stream_state.solw); + stream_state.solv, stream_state.solw, + stream_state.mesh_quad.get()); } } } diff --git a/lib/stream_reader.cpp b/lib/stream_reader.cpp index 8bc8c475..92ecb32b 100644 --- a/lib/stream_reader.cpp +++ b/lib/stream_reader.cpp @@ -743,7 +743,7 @@ void StreamState::ResetMeshAndSolution(VisualizationScene* vs) VisualizationSceneVector3d *vss = dynamic_cast(vs); - vss->NewMeshAndSolution(mesh.get(), grid_f.get()); + vss->NewMeshAndSolution(mesh.get(), mesh_quad.get(), grid_f.get()); } } } diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index 623a7497..244ed611 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -340,9 +340,10 @@ void VisualizationSceneVector3d::ToggleScalarFunction() } VisualizationSceneVector3d::VisualizationSceneVector3d(Mesh &m, Vector &sx, - Vector &sy, Vector &sz) + Vector &sy, Vector &sz, Mesh *mc) { mesh = &m; + mesh_coarse = mc; solx = &sx; soly = &sy; solz = &sz; @@ -355,7 +356,8 @@ VisualizationSceneVector3d::VisualizationSceneVector3d(Mesh &m, Vector &sx, Init(); } -VisualizationSceneVector3d::VisualizationSceneVector3d(GridFunction &vgf) +VisualizationSceneVector3d::VisualizationSceneVector3d(GridFunction &vgf, + Mesh *mc) { FiniteElementSpace *fes = vgf.FESpace(); if (fes == NULL || fes->GetVDim() != 3) @@ -367,6 +369,7 @@ VisualizationSceneVector3d::VisualizationSceneVector3d(GridFunction &vgf) VecGridF = &vgf; mesh = fes->GetMesh(); + mesh_coarse = mc; sfes = new FiniteElementSpace(mesh, fes->FEColl(), 1, fes->GetOrdering()); GridF = new GridFunction(sfes); @@ -452,7 +455,7 @@ VisualizationSceneVector3d::~VisualizationSceneVector3d() } void VisualizationSceneVector3d::NewMeshAndSolution( - Mesh *new_m, GridFunction *new_v) + Mesh *new_m, Mesh *new_mc, GridFunction *new_v) { delete sol; if (VecGridF) @@ -487,6 +490,7 @@ void VisualizationSceneVector3d::NewMeshAndSolution( VecGridF = new_v; mesh = new_m; + mesh_coarse = new_mc; FindNodePos(); sfes = new FiniteElementSpace(mesh, new_fes->FEColl(), 1, @@ -980,15 +984,23 @@ void VisualizationSceneVector3d::PrepareLines() switch (drawmesh) { case 1: - line.glBegin(GL_LINE_LOOP); - - for (j = 0; j < pointmat.Size(); j++) + { + if (mesh_coarse) { - line.glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j)); + DrawBdrElCoarseSurfEdges(line, i, pointmat); } - line.glEnd(); - break; + else + { + line.glBegin(GL_LINE_LOOP); + for (j = 0; j < pointmat.Size(); j++) + { + line.glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j)); + } + line.glEnd(); + } + break; + } case 2: for (j = 0; j < pointmat.Size(); j++) { @@ -1007,7 +1019,7 @@ void VisualizationSceneVector3d::PrepareLines() void VisualizationSceneVector3d::PrepareLines2() { - int i, j, k, fn, fo, di = 0; + int fn, fo, di = 0; double bbox_diam; line_buf.clear(); @@ -1026,7 +1038,7 @@ void VisualizationSceneVector3d::PrepareLines2() (bb.z[1]-bb.z[0])*(bb.z[1]-bb.z[0]) ); double sc = FaceShiftScale * bbox_diam; - for (i = 0; i < ne; i++) + for (int i = 0; i < ne; i++) { int attr = (dim == 3) ? mesh->GetBdrAttribute(i) : mesh->GetAttribute(i); if (!bdr_attr_to_show[attr-1]) { continue; } @@ -1105,29 +1117,53 @@ void VisualizationSceneVector3d::PrepareLines2() if (drawmesh == 1) { Array &REdges = RefG->RefEdges; - line.glBegin (GL_LINES); - for (k = 0; k < REdges.Size()/2; k++) + if (mesh_coarse) { - int *RE = &(REdges[2*k]); - if (sc == 0.0) { - for (j = 0; j < 2; j++) - line.glVertex3d (pointmat(0, RE[j]), pointmat(1, RE[j]), - pointmat(2, RE[j])); + DrawBdrElCoarseSurfEdges(line, i, pointmat, &RefG->RefPts, &REdges); } else { - for (j = 0; j < 2; j++) + DenseMatrix pointmat_shift = pointmat; + for (int j = 0; j < REdges.Size(); j++) { - double val = sc * (values(RE[j]) - minv) / (maxv - minv); - line.glVertex3d (pointmat(0, RE[j])+val*norm[0], - pointmat(1, RE[j])+val*norm[1], - pointmat(2, RE[j])+val*norm[2]); + double val = sc * (values(REdges[j]) - minv) / (maxv - minv); + for (int d = 0; d < 3; d++) + { + pointmat_shift(d, REdges[j]) += val*norm[d]; + } } + DrawBdrElCoarseSurfEdges(line, i, pointmat_shift, &RefG->RefPts, &REdges); } } - line.glEnd(); + else + { + line.glBegin (GL_LINES); + for (int k = 0; k < REdges.Size()/2; k++) + { + int *RE = &(REdges[2*k]); + + if (sc == 0.0) + { + for (int j = 0; j < 2; j++) + line.glVertex3d (pointmat(0, RE[j]), + pointmat(1, RE[j]), + pointmat(2, RE[j])); + } + else + { + for (int j = 0; j < 2; j++) + { + double val = sc * (values(RE[j]) - minv) / (maxv - minv); + line.glVertex3d (pointmat(0, RE[j])+val*norm[0], + pointmat(1, RE[j])+val*norm[1], + pointmat(2, RE[j])+val*norm[2]); + } + } + } + line.glEnd(); + } } else if (drawmesh == 2) { @@ -1145,13 +1181,13 @@ void VisualizationSceneVector3d::PrepareLines2() sides = 4; break; } - for (k = 0; k < RefG->RefGeoms.Size()/sides; k++) + for (int k = 0; k < RefG->RefGeoms.Size()/sides; k++) { RG = &(RefG->RefGeoms[k*sides]); if (sc == 0.0) { - for (j = 0; j < sides; j++) + for (int j = 0; j < sides; j++) { for (int ii = 0; ii < 3; ii++) { @@ -1162,7 +1198,7 @@ void VisualizationSceneVector3d::PrepareLines2() } else { - for (j = 0; j < sides; j++) + for (int j = 0; j < sides; j++) { double val = (values(RG[j]) - minv) / (maxv - minv); val *= sc; diff --git a/lib/vsvector3d.hpp b/lib/vsvector3d.hpp index 785b9bbe..204bc735 100644 --- a/lib/vsvector3d.hpp +++ b/lib/vsvector3d.hpp @@ -36,10 +36,11 @@ class VisualizationSceneVector3d : public VisualizationSceneSolution3d public: int ianim, ianimd, ianimmax, drawdisp; - VisualizationSceneVector3d(Mesh & m, Vector & sx, Vector & sy, Vector & sz); - VisualizationSceneVector3d (GridFunction &vgf); + VisualizationSceneVector3d(Mesh & m, Vector & sx, Vector & sy, Vector & sz, + Mesh *mc = NULL); + VisualizationSceneVector3d (GridFunction &vgf, Mesh *mc = NULL); - void NewMeshAndSolution(Mesh *new_m, GridFunction *new_v); + void NewMeshAndSolution(Mesh *new_m, Mesh *new_mc, GridFunction *new_v); virtual ~VisualizationSceneVector3d(); From 95f8336f485c2ace2f07e0742cba8f2f45045b58 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 19 Jul 2024 07:29:24 -0700 Subject: [PATCH 56/57] Minor optimization. --- lib/vssolution3d.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index b045a530..bc135ab2 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -1456,8 +1456,8 @@ void VisualizationSceneSolution3d::DrawCoarseSurfEdges( } auto &ref = mesh->GetRefinementTransforms(); IsoparametricTransformation trans; - BiLinear2DFiniteElement fe_face; - TriLinear3DFiniteElement fe; + static const BiLinear2DFiniteElement fe_face; + static const TriLinear3DFiniteElement fe; trans.SetFE(&fe); DenseMatrix emb_pointmat; @@ -1522,7 +1522,7 @@ void VisualizationSceneSolution3d::DrawCoarseSurfEdges( if ((emb_ip1(d) != 0. && emb_ip1(d) != 1.) || (emb_ip2(d) != 0. && emb_ip2(d) != 1.)) { inter++; } - if (e2>= 0 && ref.embeddings[e1].parent == ref.embeddings[e2].parent) + if (e2 >= 0 && ref.embeddings[e1].parent == ref.embeddings[e2].parent) { inter--; } if (inter != 1) { continue; } From f560f044a770cc2ed1ef64af0da2f5b8753635ca Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 19 Jul 2024 07:41:46 -0700 Subject: [PATCH 57/57] Updated CHANGELOG. --- CHANGELOG | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG b/CHANGELOG index 8010c9d7..e0b0e9a5 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -31,7 +31,9 @@ Version 4.2.1 (development) options are provided: piece-wise constants on a refined mesh (LOR), L2 field with DOFs collocated (interpolation), or projection to discontinuous elements (L2 projection). Use 'Q' to switch between them. High-order quadrature data is - supported only for tensor finite elements with the first two options. + supported only for tensor finite elements with the first two options. With the + first option, only the mesh lines of the original mesh are visualized. This + feature is also supported for the element-wise cutting plane (cplane=2). Version 4.2 released on May 23, 2022