diff --git a/CHANGELOG b/CHANGELOG new file mode 100644 index 0000000..49b224f --- /dev/null +++ b/CHANGELOG @@ -0,0 +1,43 @@ + GLVis visualization tool + + _/_/_/ _/ _/ _/ _/ + _/ _/ _/ _/ _/_/_/ + _/ _/_/ _/ _/ _/ _/ _/_/ + _/ _/ _/ _/ _/ _/ _/_/ + _/_/_/ _/_/_/_/ _/ _/ _/_/_/ + + http://glvis.googlecode.com + + +Version 1.1, released on Sep 13, 2010 +===================================== + +- Anti-aliasing (key 'A') now uses the OpenGL ARB_multisample extension (when + available). By default, 4x multisampling is used. See file lib/tk.cpp. This + value can be changed by setting GLVIS_MULTISAMPLE during compilation. + +- When drawing subdivided elements, the real normals are now used (at least in + some cases) to smooth the appearance inside the element, i.e. the surface in + 2D and the (curved) boundary in 3D. + +- Enabled the shrinking of elements (F3/F4) and material subdomains (F11/F12) + for 3D meshes saved using Mesh::PrintWithPartitioning(). This allows for + better visualization of the interior of the mesh. + +- Improved the hex-cutting algorithm to handle all cases of non-flat faces. + +- Scripts now work with scalar 3D mesh/solution. + +- Changed the makefile to recompile only files that have been changed and to + allow 'make -j'. + +- Various small fixes and styling updates. + + +Version 1.0, released on Jul 21, 2010 +===================================== + +- Uploaded to http://glvis.googlecode.com. + +- Initial release. + diff --git a/INSTALL b/INSTALL index 752d05f..91d1223 100644 --- a/INSTALL +++ b/INSTALL @@ -1,5 +1,5 @@ GLVis visualization tool - version 1.0 + version 1.1 _/_/_/ _/ _/ _/ _/ _/ _/ _/ _/ _/_/_/ diff --git a/README b/README index 5dbaea8..424dcae 100644 --- a/README +++ b/README @@ -1,5 +1,5 @@ GLVis visualization tool - version 1.0 + version 1.1 _/_/_/ _/ _/ _/ _/ _/ _/ _/ _/ _/_/_/ @@ -205,5 +205,5 @@ v/V - Add/Delete a level surface value w/W - Move boundary elements up/down in direction of the normal (i.e. "plot" the boundary values in normal direction) -F3/F4 - Shrink/Zoom subdomain in a disjoint mesh -F11/F12 - Shrink/Zoom elements in a disjoint mesh +F3/F4 - Shrink/Zoom boundary elements (to the centers of the attributes) +F11/F12 - Shrink/Zoom materials subdomains (to the centers of the attributes) diff --git a/SConstruct b/SConstruct index 80c07f0..f0cec16 100644 --- a/SConstruct +++ b/SConstruct @@ -23,6 +23,9 @@ env = Environment() CC_OPTS = '-O3' DEBUG_OPTS = '-g -DGLVIS_DEBUG -Wall' +# GLVis-specific options +env.Append(CPPDEFINES = ['GLVIS_MULTISAMPLE=4']) + # Debug options debug = ARGUMENTS.get('debug', 0) if int(debug): diff --git a/glvis.cpp b/glvis.cpp index 10af17d..e991e56 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -13,9 +13,11 @@ // GLVis - an OpenGL visualization server based on the MFEM library +#include #include #include #include +#include #include "mfem.hpp" #include "lib/visual.hpp" @@ -47,6 +49,8 @@ const char *window_titles[2] ={ "GLVis [scalar data]", "GLVis [vector data]" }; istream *script = NULL; int scr_sol_autoscale = 0; +int scr_running = 0; +int scr_level = 0; Vector *init_nodes = NULL; double scr_min_val, scr_max_val; @@ -56,128 +60,165 @@ void ReadSerial(); // read the mesh and the solution from multiple files void ReadParallel(); -void InitVis(int t) +int InitVis(int t) { - InitVisualization(window_titles[t], window_x, window_y, window_w, window_h); + return InitVisualization(window_titles[t], window_x, window_y, + window_w, window_h); } // Visualize the content of a string stream (e.g. from socket/file) void VisStream(istream * iss, char * data_type) { + int window_err = 0; + if (strcmp(data_type,"fem2d_data") == 0) { - InitVis(0); - mesh = new Mesh(*iss, 0, 0); - sol.Load(*iss, mesh->GetNV()); + window_err = InitVis(0); + if (!window_err) + { + mesh = new Mesh(*iss, 0, 0); + sol.Load(*iss, mesh->GetNV()); - vs = new VisualizationSceneSolution(*mesh, sol); + vs = new VisualizationSceneSolution(*mesh, sol); + } } else if (strcmp(data_type,"vfem2d_data") == 0 || strcmp(data_type,"vfem2d_data_keys") == 0 ) { - InitVis(1); - mesh = new Mesh(*iss, 0, 0); - solu.Load(*iss, mesh->GetNV()); - solv.Load(*iss, mesh->GetNV()); + window_err = InitVis(1); + if (!window_err) + { + mesh = new Mesh(*iss, 0, 0); + solu.Load(*iss, mesh->GetNV()); + solv.Load(*iss, mesh->GetNV()); - vs = new VisualizationSceneVector(*mesh, solu, solv); + vs = new VisualizationSceneVector(*mesh, solu, solv); - if (strcmp(data_type,"vfem2d_data_keys") == 0) - *iss >> keys; + if (strcmp(data_type,"vfem2d_data_keys") == 0) + *iss >> keys; + } } else if (strcmp(data_type,"fem3d_data") == 0) { - InitVis(0); - mesh = new Mesh(*iss, 0, 0); - sol.Load(*iss, mesh->GetNV()); + window_err = InitVis(0); + if (!window_err) + { + mesh = new Mesh(*iss, 0, 0); + sol.Load(*iss, mesh->GetNV()); - vs = new VisualizationSceneSolution3d(*mesh, sol); + vs = new VisualizationSceneSolution3d(*mesh, sol); + } } else if (strcmp(data_type,"vfem3d_data") == 0 || strcmp(data_type,"vfem3d_data_keys") == 0 ) { - InitVis(1); - mesh = new Mesh(*iss, 0, 0); - solu.Load(*iss, mesh->GetNV()); - solv.Load(*iss, mesh->GetNV()); - solw.Load(*iss, mesh->GetNV()); + window_err = InitVis(1); + if (!window_err) + { + mesh = new Mesh(*iss, 0, 0); + solu.Load(*iss, mesh->GetNV()); + solv.Load(*iss, mesh->GetNV()); + solw.Load(*iss, mesh->GetNV()); - vs = new VisualizationSceneVector3d(*mesh, solu, solv, solw); + vs = new VisualizationSceneVector3d(*mesh, solu, solv, solw); - if (strcmp(data_type,"vfem3d_data_keys") == 0) - *iss >> keys; + if (strcmp(data_type,"vfem3d_data_keys") == 0) + *iss >> keys; + } } else if (strcmp(data_type,"fem3d_paragrid") == 0) { - InitVis(0); - vs = new VisualizationSceneSolution3d(*mesh, *data[0]); + window_err = InitVis(0); + if (!window_err) + vs = new VisualizationSceneSolution3d(*mesh, *data[0]); } else if (strcmp(data_type,"vfem3d_paragrid") == 0) { - InitVis(1); - vs = new VisualizationSceneVector3d(*mesh, *data[0], - *data[1], *data[2]); + window_err = InitVis(1); + if (!window_err) + vs = new VisualizationSceneVector3d(*mesh, *data[0], + *data[1], *data[2]); } else if (strcmp(data_type,"fem2d_gf_data") == 0 || strcmp(data_type,"fem2d_gf_data_keys") == 0) { - InitVis(0); - mesh = new Mesh(*iss, 1, 0); - *iss >> ws; - grid_f = new GridFunction (mesh, *iss); - grid_f -> GetNodalValues (sol); - VisualizationSceneSolution *vss; - vs = vss = new VisualizationSceneSolution(*mesh, sol); - vss -> SetGridFunction (*grid_f); - if (strcmp(data_type,"fem2d_gf_data_keys") == 0) - *iss >> keys; + window_err = InitVis(0); + if (!window_err) + { + mesh = new Mesh(*iss, 1, 0); + *iss >> ws; + grid_f = new GridFunction (mesh, *iss); + grid_f -> GetNodalValues (sol); + VisualizationSceneSolution *vss; + vs = vss = new VisualizationSceneSolution(*mesh, sol); + vss -> SetGridFunction (*grid_f); + if (strcmp(data_type,"fem2d_gf_data_keys") == 0) + *iss >> keys; + } } else if (strcmp(data_type,"vfem2d_gf_data") == 0 || strcmp(data_type,"vfem2d_gf_data_keys") == 0 ) { - InitVis(1); - mesh = new Mesh(*iss, 1, 0); - *iss >> ws; - grid_f = new GridFunction (mesh, *iss); - vs = new VisualizationSceneVector (*grid_f); - if (strcmp(data_type,"vfem2d_gf_data_keys") == 0) - *iss >> keys; + window_err = InitVis(1); + if (!window_err) + { + mesh = new Mesh(*iss, 1, 0); + *iss >> ws; + grid_f = new GridFunction (mesh, *iss); + vs = new VisualizationSceneVector (*grid_f); + if (strcmp(data_type,"vfem2d_gf_data_keys") == 0) + *iss >> keys; + } } else if (strcmp(data_type,"fem3d_gf_data") == 0 || strcmp(data_type,"fem3d_gf_data_keys") == 0) { - InitVis(0); - mesh = new Mesh(*iss, 1, 0); - *iss >> ws; - grid_f = new GridFunction (mesh, *iss); - grid_f -> GetNodalValues (sol); - VisualizationSceneSolution3d *vss; - vs = vss = new VisualizationSceneSolution3d(*mesh, sol); - vss -> SetGridFunction (grid_f); - if (strcmp(data_type,"fem3d_gf_data_keys") == 0) - *iss >> keys; + window_err = InitVis(0); + if (!window_err) + { + mesh = new Mesh(*iss, 1, 0); + *iss >> ws; + grid_f = new GridFunction (mesh, *iss); + grid_f -> GetNodalValues (sol); + VisualizationSceneSolution3d *vss; + vs = vss = new VisualizationSceneSolution3d(*mesh, sol); + vss -> SetGridFunction (grid_f); + if (strcmp(data_type,"fem3d_gf_data_keys") == 0) + *iss >> keys; + } } else if (strcmp(data_type,"vfem3d_gf_data") == 0 || - strcmp(data_type,"vfem3d_gf_data_keys") == 0 ) + strcmp(data_type,"vfem3d_gf_data_keys") == 0) { - InitVis(1); - mesh = new Mesh(*iss, 1, 0); - *iss >> ws; - grid_f = new GridFunction (mesh, *iss); - vs = new VisualizationSceneVector3d (*grid_f); - if (strcmp(data_type,"vfem3d_gf_data_keys") == 0) - *iss >> keys; + window_err = InitVis(1); + if (!window_err) + { + mesh = new Mesh(*iss, 1, 0); + *iss >> ws; + grid_f = new GridFunction (mesh, *iss); + vs = new VisualizationSceneVector3d (*grid_f); + if (strcmp(data_type,"vfem3d_gf_data_keys") == 0) + *iss >> keys; + } } - else{ + else + { cerr << "Unknown data format" << endl; cerr << data_type << endl; return; } - SetVisualizationScene (vs, 3, keys); - KillVisualization(); // deletes vs - vs = NULL; - delete grid_f; grid_f = NULL; - delete mesh; mesh = NULL; + if (!window_err) + { + SetVisualizationScene (vs, 3, keys); + KillVisualization(); // deletes vs + vs = NULL; + delete grid_f; grid_f = NULL; + delete mesh; mesh = NULL; + } + else + { + cerr << "Initializing the visualization failed." << endl; + } } int ScriptReadSolution(istream &scr, Mesh **mp, GridFunction **sp) @@ -289,13 +330,14 @@ void ExecuteScriptCommand() #endif string word; - int done_one_command = 0, running = 0; - while (running || !done_one_command) + int done_one_command = 0; + while (!done_one_command) { scr >> ws; if (!scr.good()) { cout << "End of script." << endl; + scr_level = 0; return; } if (scr.peek() == '#') @@ -306,12 +348,13 @@ void ExecuteScriptCommand() scr >> word; if (word == "{") { - running++; - continue; + scr_level++; } else if (word == "}") { - running--; + scr_level--; + if (scr_level < 0) + scr_level = 0; } else if (word == "solution" || word == "mesh") { @@ -363,7 +406,18 @@ void ExecuteScriptCommand() } else { - // 3D ... + if (new_g->VectorDim() == 1) + { + VisualizationSceneSolution3d *vss = + dynamic_cast(vs); + new_g->GetNodalValues(sol); + vss->NewMeshAndSolution(new_m, &sol, new_g, + scr_sol_autoscale); + } + else + { + // 3D vector field ... + } } delete grid_f; grid_f = new_g; delete mesh; mesh = new_m; @@ -547,6 +601,15 @@ void ExecuteScriptCommand() SendKeySequence(keys); MyExpose(); } + else if (word == "palette") + { + int pal; + scr >> pal; + cout << "Script: palette: " << pal << endl; + Set_Palette(pal); + vs->EventUpdateColors(); + MyExpose(); + } else { cout << "Unknown command in script: " << word << endl; @@ -556,6 +619,29 @@ void ExecuteScriptCommand() } } +void ScriptControl(); + +void ScriptIdleFunc() +{ + ExecuteScriptCommand(); + if (scr_level == 0) + ScriptControl(); +} + +void ScriptControl() +{ + if (scr_running) + { + scr_running = 0; + RemoveIdleFunc(ScriptIdleFunc); + } + else + { + scr_running = 1; + AddIdleFunc(ScriptIdleFunc); + } +} + void PlayScript(istream &scr) { string word; @@ -603,10 +689,16 @@ void PlayScript(istream &scr) } } + int window_err; if (grid_f->VectorDim() == 1) - InitVis(0); + window_err = InitVis(0); else - InitVis(1); + window_err = InitVis(1); + if (window_err) + { + cerr << "Initializing the visualization failed." << endl; + return; + } if (mesh->Dimension() == 2) { @@ -624,10 +716,23 @@ void PlayScript(istream &scr) } else if (mesh->Dimension() == 3) { + if (grid_f->VectorDim() == 1) + { + grid_f->GetNodalValues(sol); + VisualizationSceneSolution3d *vss; + vs = vss = new VisualizationSceneSolution3d(*mesh, sol); + vss->SetGridFunction(grid_f); + } + else + { + // ... + } } - auxKeyFunc(XK_space, ExecuteScriptCommand); + scr_level = scr_running = 0; + auxKeyFunc(XK_space, ScriptControl); script = &scr; + scr_level = scr_running = 0; SetVisualizationScene (vs, 3); KillVisualization(); // deletes vs @@ -708,11 +813,11 @@ int main (int argc, char *argv[]) } cout << endl - << " _/_/_/ _/ _/ _/ _/ " << endl - << " _/ _/ _/ _/ _/_/_/ " << endl - << " _/ _/_/ _/ _/ _/ _/ _/_/ " << endl - << " _/ _/ _/ _/ _/ _/ _/_/ " << endl - << " _/_/_/ _/_/_/_/ _/ _/ _/_/_/ " << endl + << " _/_/_/ _/ _/ _/ _/" << endl + << " _/ _/ _/ _/ _/_/_/" << endl + << " _/ _/_/ _/ _/ _/ _/ _/_/" << endl + << " _/ _/ _/ _/ _/ _/ _/_/" << endl + << " _/_/_/ _/_/_/_/ _/ _/ _/_/_/" << endl << endl ; // get rid of zombies @@ -879,43 +984,50 @@ int main (int argc, char *argv[]) else ReadSerial(); + int window_err; if (mesh->Dimension() == 2) { if ((input & 8) == 0) { VisualizationSceneSolution *vss; - InitVis(0); - if ((input & 4) == 0) - Set_Palette(4); // Set_Palette(11); - vs = vss = new VisualizationSceneSolution(*mesh, sol); - if (is_gf) - { - vss->SetGridFunction(*grid_f); - vss->ToggleShading(); - } - if ((input & 4) == 0) + window_err = InitVis(0); + if (!window_err) { - vs->OrthogonalProjection = 1; - vs->light = 0; - vs->Zoom(1.8); - // increase the refinement factors if the mesh is "curved" - if (mesh->GetNodalFESpace()) + if ((input & 4) == 0) + Set_Palette(4); // Set_Palette(11); + vs = vss = new VisualizationSceneSolution(*mesh, sol); + if (is_gf) + { + vss->SetGridFunction(*grid_f); + vss->ToggleShading(); + } + if ((input & 4) == 0) { - cout << "Curvilinear mesh, subdivision factors = 4, 1" - << endl; - vs->SetRefineFactors(4, 1); + vs->OrthogonalProjection = 1; + vs->light = 0; + vs->Zoom(1.8); + // increase the refinement factors if the mesh is "curved" + if (mesh->GetNodalFESpace()) + { + cout << "Curvilinear mesh, subdivision factors = 4, 1" + << endl; + vs->SetRefineFactors(4, 1); + } + vs->SetValueRange(-1.0*(grid_f->Max()+1.0), + 1.0*(grid_f->Max()+1.0)); } - vs->SetValueRange(-1.0*(grid_f->Max()+1.0), - 1.0*(grid_f->Max()+1.0)); } } else { - InitVis(1); - if (is_gf) - vs = new VisualizationSceneVector (*grid_f); - else - vs = new VisualizationSceneVector(*mesh, solu, solv); + window_err = InitVis(1); + if (!window_err) + { + if (is_gf) + vs = new VisualizationSceneVector(*grid_f); + else + vs = new VisualizationSceneVector(*mesh, solu, solv); + } } } else // 3D @@ -923,42 +1035,56 @@ int main (int argc, char *argv[]) if ((input & 8) == 0 && (input & 512) == 0) { VisualizationSceneSolution3d *vss; - InitVis(0); - vs = vss = new VisualizationSceneSolution3d(*mesh, sol); - if (is_gf) - { - vss->SetGridFunction (grid_f); - vss->ToggleShading(); - vss->Prepare(); - // vss -> Draw(); // ??? - } - if ((input & 4) == 0) + window_err = InitVis(0); + if (!window_err) { - // Set_Palette(4); - Set_Palette(11); - Set_Material_And_Light(4,3); - vss->ToggleDrawAxes(); - vss->ToggleDrawMesh(); - vs->SetValueRange(-1.0*(grid_f->Max()+1.0), - 1.0*(grid_f->Max()+1.0)); + vs = vss = new VisualizationSceneSolution3d(*mesh, sol); + if (is_gf) + { + vss->SetGridFunction (grid_f); + vss->ToggleShading(); + vss->Prepare(); + // vss -> Draw(); // ??? + } + if ((input & 4) == 0) + { + // Set_Palette(4); + Set_Palette(11); + Set_Material_And_Light(4,3); + vss->ToggleDrawAxes(); + vss->ToggleDrawMesh(); + vs->SetValueRange(-1.0*(grid_f->Max()+1.0), + 1.0*(grid_f->Max()+1.0)); + } } } else { - InitVis(1); - if (is_gf) - vs = new VisualizationSceneVector (*grid_f); - else - vs = new VisualizationSceneVector3d(*mesh, solu, solv, solw); + window_err = InitVis(1); + if (!window_err) + { + if (is_gf) + vs = new VisualizationSceneVector3d(*grid_f); + else + vs = new VisualizationSceneVector3d(*mesh, solu, solv, solw); + } } } - if (mesh->Dimension() == 2 && (input & 12) == 0) - SetVisualizationScene(vs, 2, keys); + if (!window_err) + { + if (mesh->Dimension() == 2 && (input & 12) == 0) + SetVisualizationScene(vs, 2, keys); + else + SetVisualizationScene(vs, 3, keys); + KillVisualization(); // deletes vs + if (is_gf) delete grid_f; + delete mesh; + } else - SetVisualizationScene(vs, 3, keys); - KillVisualization(); // deletes vs - if (is_gf) delete grid_f; - delete mesh; + { + cerr << "Initializing the visualization failed." << endl; + return 1; + } } cout << "Thank you for using GLVis." << endl; @@ -1071,10 +1197,12 @@ void ReadParallel() istream ** parin = new istream* [np]; int * dim = new int[np]; - for (i = 0; i < np; i++) { + for (i = 0; i < np; i++) + { sprintf(fname,"%s_%d",mesh_file,i); parin[i] = new ifstream(fname); - if (!(*parin[i])) { + if (!(*parin[i])) + { cerr << "Can not open mesh file " << fname << ". Exit. \n"; exit(1); } @@ -1087,22 +1215,26 @@ void ReadParallel() meshout.precision(12); mesh -> Print(meshout); - for (i = 0; i < np; i++) { + for (i = 0; i < np; i++) + { sprintf(fname,"%s_%d",sol_file,i); parin[i] = new ifstream(fname); - if (!(*parin[i])) { + if (!(*parin[i])) + { cerr << "Can not open solution file " << fname << ". Exit. \n"; exit(1); } } - if (input & 128) { + if (input & 128) + { // get rid of NetGen's info line for (i = 0; i < np; i++) parin[i] -> getline(fname,128); sol.Load(parin,np,dim); } - else if (input & 512) { + else if (input & 512) + { solu.Load(parin, np, dim); solv.Load(parin, np, dim); solw.Load(parin, np, dim); diff --git a/lib/aux_gl.cpp b/lib/aux_gl.cpp index 9aeb1f4..6c1bd80 100644 --- a/lib/aux_gl.cpp +++ b/lib/aux_gl.cpp @@ -289,10 +289,13 @@ GLenum auxInitWindow(const char *title) { int useDoubleAsSingle = 0; - if (tkInitWindow(title) == GL_FALSE) { - if (AUX_WIND_IS_SINGLE(displayModeType)) { + if (tkInitWindow(title) == GL_FALSE) + { + if (AUX_WIND_IS_SINGLE(displayModeType)) + { tkInitDisplayMode(displayModeType | AUX_DOUBLE); - if (tkInitWindow(title) == GL_FALSE) { + if (tkInitWindow(title) == GL_FALSE) + { return GL_FALSE; /* curses, foiled again */ } fprintf(stderr, "Can't initialize a single buffer visual.\n"); @@ -301,6 +304,10 @@ GLenum auxInitWindow(const char *title) displayModeType = displayModeType | AUX_DOUBLE; useDoubleAsSingle = 1; } + else + { + return GL_FALSE; + } } tkReshapeFunc(DefaultHandleReshape); tkExposeFunc(DefaultHandleExpose); diff --git a/lib/aux_vis.cpp b/lib/aux_vis.cpp index 03317bf..e0ab319 100644 --- a/lib/aux_vis.cpp +++ b/lib/aux_vis.cpp @@ -13,6 +13,8 @@ #include #include +#include "mfem.hpp" + #include "palettes.hpp" #include "aux_vis.hpp" #include "gl2ps.h" @@ -34,9 +36,8 @@ float MatAlphaCenter = 0.5; void MyExpose(GLsizei w, GLsizei h); -void InitVisualization (const char name[], int x, int y, int w, int h) +int InitVisualization (const char name[], int x, int y, int w, int h) { - static int init = 0; if (!init) @@ -49,9 +50,12 @@ void InitVisualization (const char name[], int x, int y, int w, int h) cout << "OpenGL Visualization" << endl; #endif - auxInitDisplayMode (GLenum(AUX_DOUBLE | AUX_RGBA | AUX_DEPTH )); - auxInitPosition (x, y, w, h); - auxInitWindow (name); + GLenum mode = AUX_DOUBLE | AUX_RGBA | AUX_DEPTH; + // mode |= (AUX_ALPHA | AUX_ACCUM); + auxInitDisplayMode(mode); + auxInitPosition(x, y, w, h); + if (auxInitWindow(name) == GL_FALSE) + return 1; Set_Texture_Image(); @@ -147,6 +151,8 @@ void InitVisualization (const char name[], int x, int y, int w, int h) auxKeyFunc(XK_parenright, EnlargeWindow); locscene = NULL; + + return 0; } void SendKeyEvent (KeySym keysym, int Shift=0) @@ -158,10 +164,13 @@ void SendKeyEvent (KeySym keysym, int Shift=0) xke.type = KeyPress; xke.state = 0; - if (!Shift) { + if (!Shift) + { xke.keycode = XKeysymToKeycode(xke.display, keysym); XSendEvent(auxXDisplay(), auxXWindow(), True, KeyPressMask, (XEvent*)&xke); - } else { + } + else + { xke.keycode = XKeysymToKeycode(xke.display, XK_Shift_L); XSendEvent(auxXDisplay(), auxXWindow(), True, KeyPressMask, (XEvent*)&xke); xke.state |= ShiftMask; @@ -181,24 +190,28 @@ void SendKeySequence (char * seq) { char * key = seq; - for ( ; *key != '\0'; key++ ) { // see /usr/include/X11/keysymdef.h - - if ( ((*key - '0') < 10) && ((*key - '0') >= 0) ) { // (keypad) number + for ( ; *key != '\0'; key++ ) // see /usr/include/X11/keysymdef.h + { + if ( ((*key - '0') < 10) && ((*key - '0') >= 0) ) // (keypad) number + { SendKeyEvent(XK_0 + (*key) -'0'); continue; } - if ( ((*key - 'a') < 26) && ((*key - 'a') >= 0) ) { // lowercase letter + if ( ((*key - 'a') < 26) && ((*key - 'a') >= 0) ) // lowercase letter + { SendKeyEvent(XK_a + (*key) -'a'); continue; } - if ( ((*key - 'A') < 26) && ((*key - 'A') >= 0) ) { // uppercase letter + if ( ((*key - 'A') < 26) && ((*key - 'A') >= 0) ) // uppercase letter + { SendKeyEvent(XK_A + (*key) -'A',1); continue; } - switch (*key) { + switch (*key) + { case '+': SendKeyEvent(XK_plus); continue; @@ -231,7 +244,8 @@ void SendKeySequence (char * seq) continue; case '~': // special codes key++; - switch (*key) { + switch (*key) + { case 'e': // expose event SendExposeEvent(); break; @@ -271,7 +285,10 @@ void SendKeySequence (char * seq) } } -void SetVisualizationScene (VisualizationScene * scene, int view, char * keys){ +void InitIdleFuncs(); + +void SetVisualizationScene(VisualizationScene * scene, int view, char * keys) +{ locscene = scene; locscene -> view = view; @@ -280,16 +297,16 @@ void SetVisualizationScene (VisualizationScene * scene, int view, char * keys){ else scene -> CenterObject(); + InitIdleFuncs(); if (scene -> spinning) - auxIdleFunc(MainLoop); - else - auxIdleFunc((void (*)())0); + AddIdleFunc(MainLoop); if (keys) SendKeySequence(keys); - // auxMainLoop(MainLoop); auxMainLoop(NULL); + + InitIdleFuncs(); } void KillVisualization() @@ -303,7 +320,8 @@ void KillVisualization() auxCloseWindow(); } -void SendExposeEvent(){ +void SendExposeEvent() +{ XExposeEvent ev; ev.type = Expose; ev.count = 0; @@ -324,7 +342,7 @@ void MyReshape(GLsizei w, GLsizei h) double ViewCenterX = locscene->ViewCenterX; double ViewCenterY = locscene->ViewCenterY; - if ( locscene->OrthogonalProjection ) + if (locscene->OrthogonalProjection) { double scale = locscene->ViewScale; if (w <= h) @@ -351,7 +369,7 @@ void MyReshape(GLsizei w, GLsizei h) void MyExpose(GLsizei w, GLsizei h) { MyReshape (w, h); - locscene -> Draw (); + locscene -> Draw(); } void MyExpose() @@ -362,6 +380,38 @@ void MyExpose() MyExpose(wa.width, wa.height); } + +Array IdleFuncs; +int LastIdleFunc; + +void InitIdleFuncs() +{ + IdleFuncs.SetSize(0); + LastIdleFunc = 0; + auxIdleFunc(NULL); +} + +void MainIdleFunc() +{ + LastIdleFunc = (LastIdleFunc + 1) % IdleFuncs.Size(); + if (IdleFuncs[LastIdleFunc]) + (*IdleFuncs[LastIdleFunc])(); +} + +void AddIdleFunc(void (*Func)(void)) +{ + IdleFuncs.Union(Func); + auxIdleFunc(MainIdleFunc); +} + +void RemoveIdleFunc(void (*Func)(void)) +{ + IdleFuncs.DeleteFirst(Func); + if (IdleFuncs.Size() == 0) + auxIdleFunc(NULL); +} + + double xang, yang; double srot[16], sph_t, sph_u; static GLint oldx, oldy, startx, starty; @@ -369,14 +419,16 @@ static GLint oldx, oldy, startx, starty; int constrained_spinning = 0; -void MainLoop(){ +void MainLoop() +{ static int p = 1; struct timespec req; - if (locscene -> spinning){ + if (locscene -> spinning) + { if (!constrained_spinning) { glMatrixMode (GL_MODELVIEW); - glLoadIdentity (); + glLoadIdentity(); glRotatef(xang, 0.0f, 1.0f, 0.0f); glRotatef(yang, 1.0f, 0.0f, 0.0f); glMultMatrixd (locscene -> rotmat); @@ -395,7 +447,8 @@ void MainLoop(){ req.tv_nsec = 10000000; nanosleep (&req, NULL); // sleep for 0.01 seconds } - if (locscene -> movie) { + if (locscene -> movie) + { char fname[20]; sprintf(fname, "GLVis_m%04d", p++); Screenshot(fname); @@ -404,12 +457,14 @@ void MainLoop(){ //== PRESSED MOUSE BUTTONS EVENTS ============================================ -inline double sqr(double t){ +inline double sqr(double t) +{ return t*t; } inline void ComputeSphereAngles(int &newx, int &newy, - double &new_sph_u, double &new_sph_t){ + double &new_sph_u, double &new_sph_t) +{ GLint viewport[4]; double r, x, y, rr; const double maxr = 0.996194698091745532295; @@ -428,9 +483,10 @@ inline void ComputeSphereAngles(int &newx, int &newy, new_sph_t = atan2(y, x); } -void LeftButtonDown (AUX_EVENTREC *event){ +void LeftButtonDown (AUX_EVENTREC *event) +{ locscene -> spinning = 0; - auxIdleFunc((void (*)())0); + RemoveIdleFunc(MainLoop); oldx = event->data[AUX_MOUSEX]; oldy = event->data[AUX_MOUSEY]; @@ -444,7 +500,8 @@ void LeftButtonDown (AUX_EVENTREC *event){ starty = oldy; } -void LeftButtonLoc (AUX_EVENTREC *event){ +void LeftButtonLoc (AUX_EVENTREC *event) +{ GLint newx = event->data[AUX_MOUSEX]; GLint newy = event->data[AUX_MOUSEY]; int sendexpose = 1; @@ -514,16 +571,18 @@ void LeftButtonLoc (AUX_EVENTREC *event){ SendExposeEvent(); } -void LeftButtonUp (AUX_EVENTREC *event){ +void LeftButtonUp (AUX_EVENTREC *event) +{ GLint newx = event->data[AUX_MOUSEX]; GLint newy = event->data[AUX_MOUSEY]; xang = (newx-startx)/5.0; yang = (newy-starty)/5.0; - if ( (event->data[2] & ShiftMask) && (xang != 0.0 || yang != 0.0) ){ + if ( (event->data[2] & ShiftMask) && (xang != 0.0 || yang != 0.0) ) + { locscene -> spinning = 1; - auxIdleFunc(MainLoop); + AddIdleFunc(MainLoop); if (xang > 20) xang = 20; if (xang < -20) xang = -20; if (yang > 20) yang = 20; if (yang < -20) yang = -20; @@ -534,12 +593,14 @@ void LeftButtonUp (AUX_EVENTREC *event){ } } -void MiddleButtonDown (AUX_EVENTREC *event){ +void MiddleButtonDown (AUX_EVENTREC *event) +{ startx = oldx = event->data[AUX_MOUSEX]; starty = oldy = event->data[AUX_MOUSEY]; } -void MiddleButtonLoc (AUX_EVENTREC *event){ +void MiddleButtonLoc (AUX_EVENTREC *event) +{ GLint newx = event->data[AUX_MOUSEX]; GLint newy = event->data[AUX_MOUSEY]; @@ -572,26 +633,21 @@ void MiddleButtonLoc (AUX_EVENTREC *event){ } void MiddleButtonUp (AUX_EVENTREC *event) -{ -/* - GLint newx = event->data[AUX_MOUSEX]; - GLint newy = event->data[AUX_MOUSEY]; - - locscene -> Translate ((double)(newx-oldx)/200,(double)(newy-oldy)/200); - SendExposeEvent(); -*/ -} +{} -void RightButtonDown (AUX_EVENTREC *event){ +void RightButtonDown (AUX_EVENTREC *event) +{ startx = oldx = event->data[AUX_MOUSEX]; starty = oldy = event->data[AUX_MOUSEY]; } -void RightButtonLoc (AUX_EVENTREC *event){ +void RightButtonLoc (AUX_EVENTREC *event) +{ GLint newx = event->data[AUX_MOUSEX]; GLint newy = event->data[AUX_MOUSEY]; - if (event->data[2] & ShiftMask) { + if (event->data[2] & ShiftMask) + { glLoadIdentity(); // GLfloat light[] = {newx,-newy, sqrt((float)(newx*newx+newy*newy)), 0.0 }; newx -= startx; @@ -629,23 +685,28 @@ void RightButtonLoc (AUX_EVENTREC *event){ } void RightButtonUp (AUX_EVENTREC *event) -{ -/* - GLint newx = event->data[AUX_MOUSEX]; - GLint newy = event->data[AUX_MOUSEY]; - - if (!(event->data[2] & (ShiftMask|ControlMask))) - locscene -> Scale ( exp ( double (oldy-newy)/50 ) ); - SendExposeEvent(); -*/ -} +{} int Screenshot(const char *fname) { +#ifdef GLVIS_DEBUG + cout << "Screenshot: glXWaitX() ... " << flush; +#endif + glXWaitX(); +#ifdef GLVIS_DEBUG + cout << "done." << endl; +#endif + +#ifdef GLVIS_DEBUG + cout << "Screenshot: glFinish() ... " << flush; +#endif glFinish(); // events in the X event queue may not be complete // in particular ExposeEvents generated by SendExposeEvent() // or keys sent by SendKeySequence +#ifdef GLVIS_DEBUG + cout << "done." << endl; +#endif #ifdef GLVIS_USE_LIBTIFF // Save a TIFF image. This requires the libtiff library, see www.libtiff.org @@ -661,7 +722,7 @@ int Screenshot(const char *fname) XGetWindowAttributes(auxXDisplay(), auxXWindow(), &wa); int w = wa.width; int h = wa.height; - MyExpose(w,h); + // MyExpose(w,h); glReadBuffer(GL_FRONT); unsigned char** pixels = new unsigned char*[ h ]; @@ -704,18 +765,21 @@ int Screenshot(const char *fname) #endif } -void KeyS () +void KeyS() { static int p = 1; - if (locscene -> spinning) { + if (locscene -> spinning) + { locscene -> movie = 1 - locscene -> movie; if (locscene -> movie) cout << "Recording a movie (series of snapshots)..." << endl; else cout << endl; // use (ImageMagik's) convert GLVis_m* GLVis.{gif,mpg} - } else { + } + else + { cout << "Taking snapshot number " << p << "... "; char fname[20]; sprintf(fname, "GLVis_s%02d", p++); @@ -724,7 +788,7 @@ void KeyS () } } -void KeyP () +void KeyP() { int state, buffsize; FILE * fp; @@ -737,7 +801,7 @@ void KeyP () state = GL2PS_OVERFLOW; locscene -> print = 1; glGetIntegerv(GL_VIEWPORT, viewport); - while( state == GL2PS_OVERFLOW ) + while (state == GL2PS_OVERFLOW) { buffsize += 1024*1024; gl2psBeginPage ( "GLVis.eps", "GLVis", viewport, @@ -767,30 +831,58 @@ void KeyP () locscene -> Draw(); } -void KeyQPressed (){ +void KeyQPressed() +{ visualize = 0; } -void Key0Pressed(){ +void CheckSpin() +{ + if (fabs(xang) < 1.e-2) + xang = 0.; + if (xang != 0. || yang != 0.) + { + locscene->spinning = 1; + AddIdleFunc(MainLoop); + } + else + { + locscene->spinning = 0; + RemoveIdleFunc(MainLoop); + } +} + +void Key0Pressed() +{ xang--; + CheckSpin(); } -void KeyDeletePressed(){ +void KeyDeletePressed() +{ if (locscene -> spinning) + { + xang = yang = 0.; locscene -> spinning = 0; - else { - xang = 1; + RemoveIdleFunc(MainLoop); + } + else + { + xang = 1.; locscene -> spinning = 1; - auxIdleFunc(MainLoop); + AddIdleFunc(MainLoop); constrained_spinning = 1; } } -void KeyEnterPressed(){ +void KeyEnterPressed() +{ xang++; + CheckSpin(); } -void Key7Pressed(){ +void Key7Pressed() +{ glMatrixMode (GL_MODELVIEW); glLoadMatrixd (locscene -> rotmat); glRotated ( 1.0, 0.0, -1.0, 0.0 ); @@ -798,12 +890,14 @@ void Key7Pressed(){ SendExposeEvent(); } -void Key8Pressed(){ +void Key8Pressed() +{ locscene -> Rotate (0,-1); SendExposeEvent(); } -void Key9Pressed(){ +void Key9Pressed() +{ glMatrixMode (GL_MODELVIEW); glLoadMatrixd (locscene -> rotmat); glRotated ( -1.0, 1.0, 0.0, 0.0 ); @@ -811,7 +905,8 @@ void Key9Pressed(){ SendExposeEvent(); } -void Key4Pressed(){ +void Key4Pressed() +{ glMatrixMode (GL_MODELVIEW); glLoadMatrixd (locscene -> rotmat); glRotated ( -1.0, 0.0, 0.0, 1.0 ); @@ -819,7 +914,8 @@ void Key4Pressed(){ SendExposeEvent(); } -void Key5Pressed(){ +void Key5Pressed() +{ if (locscene -> view == 2) locscene -> CenterObject2D(); else @@ -827,7 +923,8 @@ void Key5Pressed(){ SendExposeEvent(); } -void Key6Pressed(){ +void Key6Pressed() +{ glMatrixMode (GL_MODELVIEW); glLoadMatrixd (locscene -> rotmat); glRotated ( 1.0, 0.0, 0.0, 1.0 ); @@ -835,7 +932,8 @@ void Key6Pressed(){ SendExposeEvent(); } -void Key1Pressed(){ +void Key1Pressed() +{ glMatrixMode (GL_MODELVIEW); glLoadMatrixd (locscene -> rotmat); glRotated ( 1.0, 1.0, 0.0, 0.0 ); @@ -843,12 +941,14 @@ void Key1Pressed(){ SendExposeEvent(); } -void Key2Pressed(){ +void Key2Pressed() +{ locscene -> Rotate (0,1); SendExposeEvent(); } -void Key3Pressed(){ +void Key3Pressed() +{ glMatrixMode (GL_MODELVIEW); glLoadMatrixd (locscene -> rotmat); glRotated ( 1.0, 0.0, 1.0, 0.0 ); @@ -856,22 +956,26 @@ void Key3Pressed(){ SendExposeEvent(); } -void KeyLeftPressed(){ +void KeyLeftPressed() +{ locscene -> Rotate (-5,0); SendExposeEvent(); } -void KeyRightPressed(){ +void KeyRightPressed() +{ locscene -> Rotate (5,0); SendExposeEvent(); } -void KeyUpPressed(){ +void KeyUpPressed() +{ locscene -> Rotate (0,-5); SendExposeEvent(); } -void KeyDownPressed(){ +void KeyDownPressed() +{ locscene -> Rotate (0,5); SendExposeEvent(); } @@ -965,7 +1069,8 @@ void MoveResizeWindow(int x, int y, int w, int h) // Draw a cone of radius 1 with base in the x-y plane and center at (0,0,2) -void Cone(){ +void Cone() +{ const int n = 8; const double step = 2*M_PI/n; const double nz = (1.0/4.0); @@ -978,7 +1083,8 @@ void Cone(){ glVertex3d(0, 0, 0); glNormal3d(1.0, 0.0, nz); glVertex3d(1, 0, -4); - for(i=1; i #include #include +#include int Num_Materials = 5; int Current_Material = 3; @@ -105,7 +106,7 @@ void Set_Material() glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, memis); } -void Set_Light () +void Set_Light() { glLoadIdentity(); @@ -295,7 +296,7 @@ void Set_Black_Material() } } -void Set_Background () +void Set_Background() { switch (Background) { @@ -308,32 +309,42 @@ void Set_Background () } } -void Toggle_Background () +void Toggle_Background() { Background = 1 - Background; } -void Set_Transparency () +void Set_Transparency() { glEnable (GL_BLEND); glDepthMask (GL_FALSE); glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); } -void Remove_Transparency () +void Remove_Transparency() { glDisable (GL_BLEND); glDepthMask (GL_TRUE); } -void Set_AntiAliasing () +void Set_AntiAliasing() { - glEnable (GL_BLEND); +#ifdef GLVIS_MULTISAMPLE + glEnable(GL_MULTISAMPLE); + +#ifdef GL_MULTISAMPLE_FILTER_HINT_NV + std::string s = (const char *)glGetString(GL_EXTENSIONS); + if (s.find("GL_NV_multisample_filter_hint") != std::string::npos) + glHint(GL_MULTISAMPLE_FILTER_HINT_NV, GL_NICEST); +#endif +#endif + + glEnable(GL_BLEND); // glDisable(GL_DEPTH_TEST); - // glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // glBlendFunc (GL_SRC_ALPHA_SATURATE, GL_ONE); - glLineWidth(1.5); + glLineWidth(1.4); glEnable(GL_POLYGON_SMOOTH); glHint(GL_POLYGON_SMOOTH_HINT, GL_NICEST); @@ -358,15 +369,26 @@ void Set_AntiAliasing () */ } -void Remove_AntiAliasing () +void Remove_AntiAliasing() { glDisable(GL_POLYGON_SMOOTH); glDisable(GL_LINE_SMOOTH); glDisable(GL_POINT_SMOOTH); glLineWidth(1.); + // glEnable(GL_DEPTH_TEST); glDisable(GL_BLEND); +#ifdef GLVIS_MULTISAMPLE +#ifdef GL_MULTISAMPLE_FILTER_HINT_NV + std::string s = (const char *)glGetString(GL_EXTENSIONS); + if (s.find("GL_NV_multisample_filter_hint") != std::string::npos) + glHint(GL_MULTISAMPLE_FILTER_HINT_NV, GL_FASTEST); +#endif + + glDisable(GL_MULTISAMPLE); +#endif + // glDisable(GL_FOG); } diff --git a/lib/material.hpp b/lib/material.hpp index 8736e78..3ee626a 100644 --- a/lib/material.hpp +++ b/lib/material.hpp @@ -11,7 +11,7 @@ extern void Set_Material(); -extern void Set_Light (); +extern void Set_Light(); extern int Next_Material_And_Light(); @@ -19,14 +19,14 @@ extern void Set_Material_And_Light(int,int); extern void Set_Black_Material(); -extern void Set_Background (); +extern void Set_Background(); -extern void Toggle_Background (); +extern void Toggle_Background(); -extern void Set_Transparency (); +extern void Set_Transparency(); -extern void Remove_Transparency (); +extern void Remove_Transparency(); -extern void Set_AntiAliasing (); +extern void Set_AntiAliasing(); -extern void Remove_AntiAliasing (); +extern void Remove_AntiAliasing(); diff --git a/lib/openglvis.cpp b/lib/openglvis.cpp index e10e10a..b04e435 100644 --- a/lib/openglvis.cpp +++ b/lib/openglvis.cpp @@ -13,10 +13,10 @@ #include #include "openglvis.hpp" -VisualizationScene :: VisualizationScene () +VisualizationScene::VisualizationScene() { glMatrixMode (GL_MODELVIEW); - glLoadIdentity (); + glLoadIdentity(); glGetDoublev (GL_MODELVIEW_MATRIX, translmat); glRotatef(-60.0, 1.0f, 0.0f, 0.0f); glRotatef(-40.0, 0.0f, 0.0f, 1.0f); @@ -30,12 +30,12 @@ VisualizationScene :: VisualizationScene () ViewCenterY = 0.0; } -VisualizationScene :: ~VisualizationScene () {;} +VisualizationScene::~VisualizationScene(){} void VisualizationScene::Rotate(double anglex, double angley) { glMatrixMode (GL_MODELVIEW); - glLoadIdentity (); + glLoadIdentity(); glRotatef(anglex, 0.0f, 1.0f, 0.0f); glRotatef(angley, 1.0f, 0.0f, 0.0f); glMultMatrixd (rotmat); @@ -45,7 +45,7 @@ void VisualizationScene::Rotate(double anglex, double angley) void VisualizationScene::Translate(double _x, double _y, double _z) { glMatrixMode (GL_MODELVIEW); - glLoadIdentity (); + glLoadIdentity(); glTranslatef (_x, -_y, _z); glMultMatrixd (translmat); glGetDoublev (GL_MODELVIEW_MATRIX, translmat); @@ -56,10 +56,11 @@ void VisualizationScene::Scale(double s) Scale (s, s, s); } -void VisualizationScene::Scale(double s1, double s2, double s3){ +void VisualizationScene::Scale(double s1, double s2, double s3) +{ /* glMatrixMode (GL_MODELVIEW); - glLoadIdentity (); + glLoadIdentity(); glScaled (s1, s2, s3); glMultMatrixd (rotmat); glGetDoublev (GL_MODELVIEW_MATRIX, rotmat); @@ -82,7 +83,7 @@ void VisualizationScene::Draw(){} void VisualizationScene::SetView(double theta, double phi) { glMatrixMode (GL_MODELVIEW); - glLoadIdentity (); + glLoadIdentity(); glGetDoublev (GL_MODELVIEW_MATRIX, translmat); glRotatef(-theta, 1.0f, 0.0f, 0.0f); diff --git a/lib/openglvis.hpp b/lib/openglvis.hpp index f5b530b..8e61038 100644 --- a/lib/openglvis.hpp +++ b/lib/openglvis.hpp @@ -23,8 +23,8 @@ class VisualizationScene double xscale, yscale, zscale; public: - VisualizationScene (); - virtual ~VisualizationScene (); + VisualizationScene(); + virtual ~VisualizationScene(); int spinning, OrthogonalProjection, print, movie; double ViewAngle, ViewScale; @@ -39,8 +39,8 @@ class VisualizationScene virtual void Prepare(); virtual void CenterObject(); virtual void CenterObject2D(); - virtual void PrepareAxes (); - virtual void Draw (); + virtual void PrepareAxes(); + virtual void Draw(); void SetView(double theta, double phi); void Zoom(double factor); diff --git a/lib/tk.cpp b/lib/tk.cpp index b96fc49..ab65512 100644 --- a/lib/tk.cpp +++ b/lib/tk.cpp @@ -34,12 +34,12 @@ * * OpenGL(TM) is a trademark of Silicon Graphics, Inc. */ + #include #include -#include -#include #include #include +#include #include "tk.h" #if defined(__cplusplus) || defined(c_plusplus) @@ -244,6 +244,7 @@ void tkCloseWindow(void) glXDestroyContext(display, context); XFreeColormap(display, colorMap); XFree((char *)visualInfo); + XCloseDisplay(display); display = 0; ExposeFunc = 0; @@ -579,31 +580,28 @@ static GLenum DoNextEvent(void) void tkExec(void) { - GLenum flag; + GLenum flag; - visualize = 1; - while (visualize) { - if (IdleFunc) { - if (IdleFunc) { - (*IdleFunc)(); - } - flag = GL_FALSE; - while (XPending(display)) { - flag |= DoNextEvent(); - } - if (flag == GL_TRUE) { - if (DisplayFunc) { - (*DisplayFunc)(); - } - } - } else { - if (DoNextEvent() == GL_TRUE) { - if (DisplayFunc) { - (*DisplayFunc)(); - } - } - } - } + visualize = 1; + while (visualize) + { + if (IdleFunc) + { + (*IdleFunc)(); + flag = GL_FALSE; + while (XPending(display)) + flag |= DoNextEvent(); + if (flag == GL_TRUE) + if (DisplayFunc) + (*DisplayFunc)(); + } + else + { + if (DoNextEvent() == GL_TRUE) + if (DisplayFunc) + (*DisplayFunc)(); + } + } } void tkExposeFunc(void (*Func)(int, int)) @@ -692,7 +690,7 @@ Window tkGetXWindow(void) static XVisualInfo *FindVisual(GLenum type) { - GLenum list[20]; + GLenum list[50]; int i; i = 0; @@ -738,6 +736,22 @@ static XVisualInfo *FindVisual(GLenum type) list[i++] = 1; } + // multisampling +#ifdef GLVIS_MULTISAMPLE +#ifdef GLX_SAMPLE_BUFFERS_ARB + std::string s = glXQueryExtensionsString(display, screen); + if (s.find("GLX_ARB_multisample") != std::string::npos) + { + list[i++] = GLX_SAMPLE_BUFFERS_ARB; + list[i++] = 1; + list[i++] = GLX_SAMPLES_ARB; + list[i++] = GLVIS_MULTISAMPLE; + } +#else +#error GLX_SAMPLE_BUFFERS_ARB is not defined! +#endif +#endif + list[i] = (int)None; return glXChooseVisual(display, screen, (int *)list); @@ -751,23 +765,49 @@ static int MakeVisualType(XVisualInfo *vi) mask = 0; glXGetConfig(display, vi, GLX_RGBA, &x); +#ifdef GLVIS_DEBUG + printf("GLX_RGBA : %d\n", x); +#endif if (x) { mask |= TK_RGB; +#ifdef GLVIS_DEBUG + glXGetConfig(display, vi, GLX_RED_SIZE, &x); + glXGetConfig(display, vi, GLX_GREEN_SIZE, &y); + glXGetConfig(display, vi, GLX_BLUE_SIZE, &z); + printf("GLX_RED_SIZE : %d\n", x); + printf("GLX_GREEN_SIZE : %d\n", y); + printf("GLX_BLUE_SIZE : %d\n", z); +#endif glXGetConfig(display, vi, GLX_ALPHA_SIZE, &x); +#ifdef GLVIS_DEBUG + printf("GLX_ALPHA_SIZE : %d\n", x); +#endif if (x > 0) { mask |= TK_ALPHA; } glXGetConfig(display, vi, GLX_ACCUM_RED_SIZE, &x); glXGetConfig(display, vi, GLX_ACCUM_GREEN_SIZE, &y); glXGetConfig(display, vi, GLX_ACCUM_BLUE_SIZE, &z); +#ifdef GLVIS_DEBUG + printf("GLX_ACCUM_RED_SIZE : %d\n", x); + printf("GLX_ACCUM_GREEN_SIZE : %d\n", y); + printf("GLX_ACCUM_BLUE_SIZE : %d\n", z); +#endif if (x > 0 && y > 0 && z > 0) { mask |= TK_ACCUM; } +#ifdef GLVIS_DEBUG + glXGetConfig(display, vi, GLX_ACCUM_ALPHA_SIZE, &x); + printf("GLX_ACCUM_ALPHA_SIZE : %d\n", x); +#endif } else { mask |= TK_INDEX; } glXGetConfig(display, vi, GLX_DOUBLEBUFFER, &x); +#ifdef GLVIS_DEBUG + printf("GLX_DOUBLEBUFFER : %d\n", x); +#endif if (x) { mask |= TK_DOUBLE; } else { @@ -775,20 +815,40 @@ static int MakeVisualType(XVisualInfo *vi) } glXGetConfig(display, vi, GLX_DEPTH_SIZE, &x); +#ifdef GLVIS_DEBUG + printf("GLX_DEPTH_SIZE : %d\n", x); +#endif if (x > 0) { mask |= TK_DEPTH; } glXGetConfig(display, vi, GLX_STENCIL_SIZE, &x); +#ifdef GLVIS_DEBUG + printf("GLX_STENCIL_SIZE : %d\n", x); +#endif if (x > 0) { mask |= TK_STENCIL; } +#ifdef GLVIS_MULTISAMPLE +#ifdef GLVIS_DEBUG +#ifdef GLX_SAMPLE_BUFFERS_ARB + glXGetConfig(display, vi, GLX_SAMPLE_BUFFERS_ARB, &x); + printf("GLX_SAMPLE_BUFFERS_ARB : %d\n", x); + glXGetConfig(display, vi, GLX_SAMPLES_ARB, &x); + printf("GLX_SAMPLES_ARB : %d\n", x); +#endif +#endif +#endif + if (glXIsDirect(display, context)) { mask |= TK_DIRECT; } else { mask |= TK_INDIRECT; } +#ifdef GLVIS_DEBUG + printf("glXIsDirect : %d\n", glXIsDirect(display, context)); +#endif return mask; } @@ -845,8 +905,8 @@ GLenum tkInitWindow(const char *title) } context = glXCreateContext(display, visualInfo, None, - (TK_IS_DIRECT(windInfo.type)) ? GL_TRUE - : GL_FALSE); + (TK_IS_DIRECT(windInfo.type)) ? GL_TRUE : + GL_FALSE); if (!context) { fprintf(stderr, "Can't create a context!\n"); return GL_FALSE; diff --git a/lib/tk.h b/lib/tk.h index c4e392d..bf3873d 100644 --- a/lib/tk.h +++ b/lib/tk.h @@ -73,8 +73,8 @@ extern "C" { #define TK_HAS_ALPHA(x) (((x) & TK_ALPHA) != 0) #define TK_HAS_DEPTH(x) (((x) & TK_DEPTH) != 0) #define TK_HAS_STENCIL(x) (((x) & TK_STENCIL) != 0) -#define TK_IS_DIRECT(x) (((x) & TK_INDIRECT) != 0) -#define TK_IS_INDIRECT(x) (((x) & TK_INDIRECT) == 0) +#define TK_IS_DIRECT(x) (((x) & TK_INDIRECT) == 0) +#define TK_IS_INDIRECT(x) (((x) & TK_INDIRECT) != 0) /* ** Event Status diff --git a/lib/vsdata.cpp b/lib/vsdata.cpp index 6ea496d..67a8018 100644 --- a/lib/vsdata.cpp +++ b/lib/vsdata.cpp @@ -28,7 +28,8 @@ using namespace std; void VisualizationSceneScalarData::Arrow3(double px, double py, double pz, double vx, double vy, double vz, double length, - double cone_scale){ + double cone_scale) +{ glPushMatrix(); double xc = 0.5*(x[0]+x[1]); @@ -111,7 +112,8 @@ void VisualizationSceneScalarData::Arrow2(double px, double py, double pz, void VisualizationSceneScalarData::Arrow(double px, double py, double pz, double vx, double vy, double vz, double length, - double cone_scale){ + double cone_scale) +{ double rhos = sqrt (vx*vx+vy*vy+vz*vz); if (rhos == 0.0) return; @@ -213,11 +215,11 @@ void VisualizationSceneScalarData::Arrow(double px, double py, double pz, glEnd(); } -void VisualizationSceneScalarData -::DrawColorBar (double minval, double maxval, - Array *level, Array *levels) +void VisualizationSceneScalarData::DrawColorBar (double minval, double maxval, + Array *level, + Array *levels) { - glMatrixMode (GL_MODELVIEW); + glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity(); @@ -249,7 +251,7 @@ ::DrawColorBar (double minval, double maxval, // glEnable (GL_COLOR_MATERIAL); glNormal3d (0, 0, 1); - Set_Material (); + Set_Material(); double eps = 1e-24; @@ -325,7 +327,6 @@ ::DrawColorBar (double minval, double maxval, double val; ostringstream * buf; - if (!level) { for (i = 0; i <= 4; i++) @@ -384,9 +385,10 @@ ::DrawColorBar (double minval, double maxval, } } - glPopAttrib (); + glPopAttrib(); - if (print) { + if (print) + { if (!level) { for (i = 0; i <= 4; i++) @@ -453,7 +455,7 @@ ::DrawColorBar (double minval, double maxval, glPopMatrix(); } -void VisualizationSceneScalarData::DrawCoordinateCross () +void VisualizationSceneScalarData::DrawCoordinateCross() { glMatrixMode (GL_PROJECTION); glPushMatrix(); @@ -500,7 +502,7 @@ void VisualizationSceneScalarData::DrawCoordinateCross () if (print) gl2psText(a_labels[2],"Times",8); glCallLists (1, GL_UNSIGNED_BYTE, a_labels[2]); - glPopAttrib (); + glPopAttrib(); glMatrixMode (GL_PROJECTION); glPopMatrix(); @@ -524,22 +526,26 @@ void VisualizationSceneScalarData::DrawCoordinateCross () VisualizationSceneScalarData * vsdata; extern VisualizationScene * locscene; -void KeyCPressed (){ +void KeyCPressed() +{ vsdata -> ToggleDrawColorbar(); SendExposeEvent(); } -void KeySPressed (){ +void KeySPressed() +{ vsdata -> ToggleScaling(); SendExposeEvent(); } -void KeyaPressed (){ +void KeyaPressed() +{ vsdata -> ToggleDrawAxes(); SendExposeEvent(); } -void KeylPressed (){ +void KeylPressed() +{ vsdata -> ToggleLight(); if (! vsdata -> light) glDisable(GL_LIGHTING); @@ -548,15 +554,17 @@ void KeylPressed (){ SendExposeEvent(); } -void KeyLPressed (){ +void KeyLPressed() +{ MySetColorLogscale = !MySetColorLogscale; - vsdata -> EventUpdateColors (); + vsdata -> EventUpdateColors(); SendExposeEvent(); } -void KeyrPressed (){ +void KeyrPressed() +{ locscene -> spinning = 0; - auxIdleFunc((void (*)())0); + RemoveIdleFunc(MainLoop); vsdata -> CenterObject(); locscene -> ViewAngle = 45.0; @@ -567,15 +575,17 @@ void KeyrPressed (){ SendExposeEvent(); } -void KeyRPressed (){ +void KeyRPressed() +{ locscene -> spinning = 0; - auxIdleFunc((void (*)())0); + RemoveIdleFunc(MainLoop); glMatrixMode (GL_MODELVIEW); - glLoadIdentity (); + glLoadIdentity(); glGetDoublev (GL_MODELVIEW_MATRIX, locscene -> translmat); Set_Light(); - switch (vsdata -> key_r_state) { + switch (vsdata -> key_r_state) + { case 0: break; @@ -613,13 +623,14 @@ void KeyRPressed (){ void Next_RGB_Palette(); int Select_New_RGB_Palette(); -void KeyPPressed (){ +void KeyPPressed() +{ Next_RGB_Palette(); - vsdata -> EventUpdateColors (); + vsdata -> EventUpdateColors(); SendExposeEvent(); } -static void KeyF5Pressed () +static void KeyF5Pressed() { int n; double min, max; @@ -652,7 +663,7 @@ void KeyF6Pressed() Select_New_RGB_Palette(); - vsdata -> EventUpdateColors (); + vsdata -> EventUpdateColors(); SendExposeEvent(); } @@ -692,18 +703,19 @@ void KeyBackslashPressed() SendExposeEvent(); } -void KeyTPressed () +void KeyTPressed() { int ml; ml = Next_Material_And_Light(); - vsdata -> EventUpdateColors (); + vsdata -> EventUpdateColors(); SendExposeEvent(); cout << "New material/light : " << ml << endl; } -void KeyGPressed (){ - Toggle_Background (); +void KeyGPressed() +{ + Toggle_Background(); SendExposeEvent(); } @@ -712,7 +724,7 @@ void KeyF1Pressed() vsdata->PrintState(); } -void KeyF2Pressed () +void KeyF2Pressed() { vsdata -> EventUpdateColors(); vsdata -> PrepareLines(); @@ -720,25 +732,26 @@ void KeyF2Pressed () SendExposeEvent(); } -void KeykPressed () +void KeykPressed() { MatAlpha -= 0.05; if (MatAlpha < 0.0) MatAlpha = 0.0; - vsdata -> EventUpdateColors (); + vsdata -> EventUpdateColors(); SendExposeEvent(); } -void KeyKPressed () +void KeyKPressed() { MatAlpha += 0.05; if (MatAlpha > 1.0) MatAlpha = 1.0; - vsdata -> EventUpdateColors (); + vsdata -> EventUpdateColors(); SendExposeEvent(); } -void KeyAPressed (){ +void KeyAPressed() +{ static int a = 0; a = 1-a; @@ -753,20 +766,20 @@ void KeyAPressed (){ } -void KeyCommaPressed () +void KeyCommaPressed() { MatAlphaCenter -= 0.25; - vsdata -> EventUpdateColors (); + vsdata -> EventUpdateColors(); SendExposeEvent(); #ifdef GLVIS_DEBUG cout << "MatAlphaCenter = " << MatAlphaCenter << endl; #endif } -void KeyLessPressed () +void KeyLessPressed() { MatAlphaCenter += 0.25; - vsdata -> EventUpdateColors (); + vsdata -> EventUpdateColors(); SendExposeEvent(); #ifdef GLVIS_DEBUG cout << "MatAlphaCenter = " << MatAlphaCenter << endl; @@ -894,15 +907,16 @@ void VisualizationSceneScalarData::ToggleTexture() } -VisualizationSceneScalarData -::VisualizationSceneScalarData(Mesh & m, Vector & s){ +VisualizationSceneScalarData::VisualizationSceneScalarData( + Mesh & m, Vector & s) +{ mesh = &m; sol = &s; Init(); } -void VisualizationSceneScalarData :: Init() +void VisualizationSceneScalarData::Init() { vsdata = this; @@ -973,7 +987,11 @@ void VisualizationSceneScalarData :: Init() glEnable(GL_AUTO_NORMAL); glEnable(GL_NORMALIZE); - // add black fog +#ifdef GLVIS_MULTISAMPLE + glDisable(GL_MULTISAMPLE); +#endif + +// add black fog // glEnable(GL_FOG); // GLfloat fogcol[4] = {0,0,0,1}; // glFogfv(GL_FOG_COLOR, fogcol); @@ -988,13 +1006,13 @@ void VisualizationSceneScalarData :: Init() PrepareAxes(); } -VisualizationSceneScalarData::~VisualizationSceneScalarData () +VisualizationSceneScalarData::~VisualizationSceneScalarData() { glDeleteLists (axeslist, 1); delete CuttingPlane; } -void VisualizationSceneScalarData::SetNewScalingFromBox () +void VisualizationSceneScalarData::SetNewScalingFromBox() { // double eps = 1e-12; double eps = 0.0; @@ -1040,19 +1058,20 @@ void VisualizationSceneScalarData::SetValueRange(double min, double max) void VisualizationSceneScalarData::ResetScaling() { - SetNewScalingFromBox (); + SetNewScalingFromBox(); } -void VisualizationSceneScalarData::CenterObject() { +void VisualizationSceneScalarData::CenterObject() +{ glMatrixMode (GL_MODELVIEW); - glLoadIdentity (); + glLoadIdentity(); glGetDoublev (GL_MODELVIEW_MATRIX, translmat); Set_Light(); // scaling = 1; - // SetNewScalingFromBox (); + // SetNewScalingFromBox(); glRotatef(-60.0, 1.0f, 0.0f, 0.0f); glRotatef(-40.0, 0.0f, 0.0f, 1.0f); @@ -1062,17 +1081,17 @@ void VisualizationSceneScalarData::CenterObject() { void VisualizationSceneScalarData::CenterObject2D() { glMatrixMode (GL_MODELVIEW); - glLoadIdentity (); + glLoadIdentity(); glGetDoublev (GL_MODELVIEW_MATRIX, translmat); Set_Light(); glGetDoublev (GL_MODELVIEW_MATRIX, rotmat); } -void VisualizationSceneScalarData::PrepareAxes(){ - +void VisualizationSceneScalarData::PrepareAxes() +{ glNewList (axeslist, GL_COMPILE); - // Set_Black_Material (); + // Set_Black_Material(); glBegin (GL_LINE_LOOP); glVertex3d(x[0], y[0], z[0]); @@ -1123,7 +1142,7 @@ void VisualizationSceneScalarData::PrepareAxes(){ glPopAttrib(); } - glEndList (); + glEndList(); } void VisualizationSceneScalarData::DrawPolygonLevelLines( @@ -1207,6 +1226,26 @@ void VisualizationSceneScalarData::PrintState() << '\n' << endl; } +void VisualizationSceneScalarData::ShrinkPoints(DenseMatrix &pointmat) +{ + int k; + double cx=0.0, cy=0.0; + + for (k = 0; k < pointmat.Width(); k++) + { + cx += pointmat(0,k); + cy += pointmat(1,k); + } + cx /= pointmat.Width(); + cy /= pointmat.Width(); + + for (k = 0; k < pointmat.Width(); k++) + { + pointmat(0,k) = shrink*pointmat(0,k) + (1-shrink)*cx; + pointmat(1,k) = shrink*pointmat(1,k) + (1-shrink)*cy; + } +} + Plane::Plane(double A,double B,double C,double D) { diff --git a/lib/vsdata.hpp b/lib/vsdata.hpp index 62cae78..d03159b 100644 --- a/lib/vsdata.hpp +++ b/lib/vsdata.hpp @@ -70,6 +70,7 @@ class VisualizationSceneScalarData : public VisualizationScene Plane * CuttingPlane; int light; int key_r_state; + double shrink; VisualizationSceneScalarData() {} VisualizationSceneScalarData (Mesh & m, Vector & s); @@ -145,6 +146,9 @@ class VisualizationSceneScalarData : public VisualizationScene void DrawRuler(); void ToggleTexture(); + + /// Shrink the set of points towards their center of gravity + void ShrinkPoints(DenseMatrix &pointmat); }; #endif diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 10a6340..ec6f745 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -20,13 +20,13 @@ #include "visual.hpp" -VisualizationSceneSolution * vssol; -extern VisualizationScene * locscene; - +VisualizationSceneSolution *vssol; +extern VisualizationScene *locscene; // Definitions of some more keys -static void SolutionKeyHPressed (){ +static void SolutionKeyHPressed() +{ cout << endl << "+------------------------------------+" << endl << "| Keys |" << endl @@ -94,7 +94,7 @@ static void SolutionKeyHPressed (){ << "+------------------------------------+" << endl; } -static void KeyF8Pressed () +static void KeyF8Pressed() { int attr; @@ -125,7 +125,7 @@ static void KeyF8Pressed () SendExposeEvent(); } -static void KeyF9Pressed () +static void KeyF9Pressed() { int attr; @@ -145,7 +145,7 @@ static void KeyF9Pressed () SendExposeEvent(); } -static void KeyF10Pressed () +static void KeyF10Pressed() { int attr; @@ -164,29 +164,33 @@ static void KeyF10Pressed () SendExposeEvent(); } -static void KeyBPressed (){ +static void KeyBPressed() +{ vssol -> ToggleDrawBdr(); SendExposeEvent(); } -static void KeyMPressed (){ +static void KeyMPressed() +{ vssol -> ToggleDrawMesh(); SendExposeEvent(); } -static void KeyEPressed (){ +static void KeyEPressed() +{ vssol -> ToggleDrawElems(); SendExposeEvent(); } -static void KeyFPressed (){ +static void KeyFPressed() +{ vssol -> ToggleShading(); SendExposeEvent(); } int refine_func = 0; -void KeyiPressed () +void KeyiPressed() { int update = 1; switch (refine_func) @@ -218,9 +222,9 @@ void KeyiPressed () } if (update && vssol -> shading == 2) { - vssol -> FindNewBox (); + vssol -> FindNewBox(); locscene -> PrepareAxes(); - vssol -> PrepareLines (); + vssol -> PrepareLines(); locscene -> Prepare(); vssol -> DefaultLevelLines(); vssol -> PrepareLevelCurves(); @@ -231,7 +235,7 @@ void KeyiPressed () << ", " << vssol -> EdgeRefineFactor << endl; } -void KeyIPressed () +void KeyIPressed() { refine_func = (refine_func+1)%4; cout << "Key 'i' will: "; @@ -308,10 +312,10 @@ static void KeyF4Pressed() } } -VisualizationSceneSolution::VisualizationSceneSolution(){} +VisualizationSceneSolution::VisualizationSceneSolution() +{} -VisualizationSceneSolution -::VisualizationSceneSolution(Mesh & m, Vector & s) +VisualizationSceneSolution::VisualizationSceneSolution(Mesh &m, Vector &s) { mesh = &m; sol = &s; @@ -340,7 +344,7 @@ void VisualizationSceneSolution::Init() el_attr_to_show[i] = 1; drawbdr = 0; - VisualizationSceneScalarData :: Init(); // Calls FindNewBox() !!! + VisualizationSceneScalarData::Init(); // Calls FindNewBox() !!! SetUseTexture(1); @@ -393,7 +397,7 @@ void VisualizationSceneSolution::Init() PrepareBoundary(); } -VisualizationSceneSolution::~VisualizationSceneSolution () +VisualizationSceneSolution::~VisualizationSceneSolution() { glDeleteLists (displlist, 1); glDeleteLists (linelist, 1); @@ -474,24 +478,42 @@ void VisualizationSceneSolution::GetRefinedValues( ShrinkPoints(tr); } -void VisualizationSceneSolution::ShrinkPoints(DenseMatrix &pointmat) +int VisualizationSceneSolution::GetRefinedValuesAndNormals( + int i, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr, + DenseMatrix &normals) { - int k; - double cx=0.0, cy=0.0; + int have_normals = 0; - for (k = 0; k < pointmat.Width(); k++) + if (drawelems < 2) { - cx += pointmat(0,k); - cy += pointmat(1,k); + rsol->GetGradients(i, ir, tr); + normals.SetSize(3, tr.Width()); + for (int j = 0; j < tr.Width(); j++) + { + normals(0, j) = -tr(0, j); + normals(1, j) = -tr(1, j); + normals(2, j) = 1.; + } + have_normals = 1; + rsol->GetValues(i, ir, vals, tr); } - cx /= pointmat.Width(); - cy /= pointmat.Width(); + else + GetRefinedDetJ(i, ir, vals, tr); - for (k = 0; k < pointmat.Width(); k++) + if (shrink != 1.0) { - pointmat(0,k) = shrink*pointmat(0,k) + (1-shrink)*cx; - pointmat(1,k) = shrink*pointmat(1,k) + (1-shrink)*cy; + ShrinkPoints(tr); + if (have_normals) + { + for (int j = 0; j < tr.Width(); j++) + { + normals(0, j) /= shrink; + normals(1, j) /= shrink; + } + } } + + return have_normals; } void VisualizationSceneSolution::SetShading(int s) @@ -698,59 +720,159 @@ void VisualizationSceneSolution::FindNewBox() SetNewScalingFromBox(); } -void VisualizationSceneSolution :: PrepareFlat (){ +void DrawTriangle(const double pts[][3], const double cv[], + const double minv, const double maxv) +{ + double nor[3]; + if (Compute3DUnitNormal(pts[0], pts[1], pts[2], nor)) + return; + glBegin(GL_TRIANGLES); + glNormal3dv(nor); + for (int j = 0; j < 3; j++) + { + MySetColor(cv[j], minv, maxv); + glVertex3dv(pts[j]); + } + glEnd(); +} + +void DrawQuad(const double pts[][3], const double cv[], + const double minv, const double maxv) +{ + double nor[3]; + if (Compute3DUnitNormal(pts[0], pts[1], pts[2], pts[3], nor)) + return; + glBegin(GL_QUADS); + glNormal3dv(nor); + for (int j = 0; j < 4; j++) + { + MySetColor(cv[j], minv, maxv); + glVertex3dv(pts[j]); + } + glEnd(); +} + +void DrawPatch(const DenseMatrix &pts, Vector &vals, DenseMatrix &normals, + const int n, const Array &ind, const double minv, + const double maxv, const int normals_opt) +{ + double na[3]; + + if (normals_opt == 1 || normals_opt == -2) + { + normals.SetSize(3, pts.Width()); + normals = 0.; + for (int i = 0; i < ind.Size(); i += n) + { + int j; + if (n == 3) + j = Compute3DUnitNormal(&pts(0, ind[i]), &pts(0, ind[i+1]), + &pts(0, ind[i+2]), na); + else + j = Compute3DUnitNormal(&pts(0, ind[i]), &pts(0, ind[i+1]), + &pts(0, ind[i+2]), &pts(0, ind[i+3]), na); + if (j == 0) + for ( ; j < n; j++) + for (int k = 0; k < 3; k++) + normals(k, ind[i+j]) += na[k]; + } + } + + if (n == 3) + glBegin(GL_TRIANGLES); + else + glBegin(GL_QUADS); + if (normals_opt != 0 && normals_opt != -1) + { + if (normals_opt > 0) + { + for (int i = 0; i < ind.Size(); i++) + { + glNormal3dv(&normals(0, ind[i])); + MySetColor(vals(ind[i]), minv, maxv); + glVertex3dv(&pts(0, ind[i])); + } + } + else + { + for (int i = ind.Size()-1; i >= 0; i--) + { + glNormal3dv(&normals(0, ind[i])); + MySetColor(vals(ind[i]), minv, maxv); + glVertex3dv(&pts(0, ind[i])); + } + } + } + else + { + for (int i = 0; i < ind.Size(); i += n) + { + int j; + if (n == 3) + j = Compute3DUnitNormal(&pts(0, ind[i]), &pts(0, ind[i+1]), + &pts(0, ind[i+2]), na); + else + j = Compute3DUnitNormal(&pts(0, ind[i]), &pts(0, ind[i+1]), + &pts(0, ind[i+2]), &pts(0, ind[i+3]), na); + if (j == 0) + { + if (normals_opt == 0) + { + glNormal3dv(na); + for ( ; j < n; j++) + { + MySetColor(vals(ind[i+j]), minv, maxv); + glVertex3dv(&pts(0, ind[i+j])); + } + } + else + { + glNormal3d(-na[0], -na[1], -na[2]); + for (j = n-1; j >= 0; j--) + { + MySetColor(vals(ind[i+j]), minv, maxv); + glVertex3dv(&pts(0, ind[i+j])); + } + } + } + } + } + glEnd(); +} + +void VisualizationSceneSolution::PrepareFlat() +{ int i, j; glNewList (displlist, GL_COMPILE); - Set_Material (); + Set_Material(); int ne = mesh -> GetNE(); DenseMatrix pointmat; Array vertices; + double pts[4][3], col[4]; for (i = 0; i < ne; i++) { if (!el_attr_to_show[mesh->GetAttribute(i)-1]) continue; - switch (mesh->GetElementType(i)) - { - case Element::TRIANGLE: - glBegin (GL_TRIANGLES); - break; - - case Element::QUADRILATERAL: - glBegin (GL_QUADS); - break; - } mesh->GetPointMatrix (i, pointmat); mesh->GetElementVertices (i, vertices); - double v10[] = { pointmat(0,1)-pointmat(0,0), - pointmat(1,1)-pointmat(1,0), - (*sol)(vertices[1])-(*sol)(vertices[0]) }; - double v21[] = { pointmat(0,2)-pointmat(0,1), - pointmat(1,2)-pointmat(1,1), - (*sol)(vertices[2])-(*sol)(vertices[1]) }; - - double norm[] = { v10[1]*v21[2]-v10[2]*v21[1], - v10[2]*v21[0]-v10[0]*v21[2], - v10[0]*v21[1]-v10[1]*v21[0] }; - double rlen = 1.0/sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]); - - glNormal3d (norm[0]*rlen, norm[1]*rlen, norm[2]*rlen); - - for (j = 0; j < pointmat.Size(); j++) + for (j = 0; j < pointmat.Width(); j++) { - MySetColor ( (*sol)(vertices[j]) , minv, maxv); - glVertex3d (pointmat.Elem(0, j), - pointmat.Elem(1, j), - (*sol)(vertices[j])); + pts[j][0] = pointmat(0, j); + pts[j][1] = pointmat(1, j); + pts[j][2] = col[j] = (*sol)(vertices[j]); } - glEnd (); + if (j == 3) + DrawTriangle(pts, col, minv, maxv); + else + DrawQuad(pts, col, minv, maxv); } - glEndList (); + glEndList(); } // determines how quads and their level lines are drawn @@ -760,65 +882,18 @@ void VisualizationSceneSolution :: PrepareFlat (){ // 2 - draw 4 triangles (split using both diagonals) const int split_quads = 0; -void DrawTriangle(const double x[], const double y[], const double z[], - const double minv, const double maxv) -{ - double u[3] = { x[1]-x[0], y[1]-y[0], z[1]-z[0] }; - double v[3] = { x[2]-x[0], y[2]-y[0], z[2]-z[0] }; - double nor[3] = { u[1]*v[2]-u[2]*v[1], - u[2]*v[0]-u[0]*v[2], - u[0]*v[1]-u[1]*v[0] }; - double l = sqrt(nor[0]*nor[0]+nor[1]*nor[1]+nor[2]*nor[2]); - - if (l > 0.0) - l = 1./l; - - glBegin(GL_TRIANGLES); - glNormal3d(nor[0]*l, nor[1]*l, nor[2]*l); - for (int j = 0; j < 3; j++) - { - MySetColor(z[j], minv, maxv); - glVertex3d(x[j], y[j], z[j]); - } - glEnd(); -} - -void DrawQuad(const double x[], const double y[], const double z[], - const double minv, const double maxv) -{ - double u[3] = { x[1]-x[0], y[1]-y[0], z[1]-z[0] }; - double v[3] = { x[2]-x[0], y[2]-y[0], z[2]-z[0] }; - double nor[3] = { u[1]*v[2]-u[2]*v[1], - u[2]*v[0]-u[0]*v[2], - u[0]*v[1]-u[1]*v[0] }; - double l = sqrt(nor[0]*nor[0]+nor[1]*nor[1]+nor[2]*nor[2]); - - if (l > 0.0) - l = 1./l; - - glBegin(GL_QUADS); - glNormal3d(nor[0]*l, nor[1]*l, nor[2]*l); - for (int j = 0; j < 4; j++) - { - MySetColor(z[j], minv, maxv); - glVertex3d(x[j], y[j], z[j]); - } - glEnd(); -} - -void VisualizationSceneSolution :: PrepareFlat2 () +void VisualizationSceneSolution::PrepareFlat2() { int i, j, k; glNewList (displlist, GL_COMPILE); - Set_Material (); + Set_Material(); int ne = mesh -> GetNE(); - DenseMatrix pointmat; + DenseMatrix pointmat, pts3d, normals; Vector values; RefinedGeometry *RefG; - double x[4], y[4], z[4]; for (i = 0; i < ne; i++) { @@ -826,10 +901,22 @@ void VisualizationSceneSolution :: PrepareFlat2 () RefG = GlobGeometryRefiner.Refine (mesh->GetElementBaseGeometry(i), TimesToRefine, EdgeRefineFactor); - GetRefinedValues (i, RefG->RefPts, values, pointmat); + j = GetRefinedValuesAndNormals(i, RefG->RefPts, values, pointmat, + normals); Array &RG = RefG->RefGeoms; int sides = mesh->GetElement(i)->GetNVertices(); +#if 1 + pts3d.SetSize(3, pointmat.Width()); + for (k = 0; k < pointmat.Width(); k++) + { + pts3d(0, k) = pointmat(0, k); + pts3d(1, k) = pointmat(1, k); + pts3d(2, k) = values(k); + } + j = (j != 0) ? 2 : 0; + DrawPatch(pts3d, values, normals, sides, RG, minv, maxv, j); +#else for (k = 0; k < RG.Size()/sides; k++) { int *ind = &RG[sides*k]; @@ -837,14 +924,14 @@ void VisualizationSceneSolution :: PrepareFlat2 () { for (j = 0; j < sides; j++) { - x[j] = pointmat(0, ind[j]); - y[j] = pointmat(1, ind[j]); - z[j] = values(ind[j]); + pts[j][0] = pointmat(0, ind[j]); + pts[j][1] = pointmat(1, ind[j]); + pts[j][2] = col[j] = values(ind[j]); } if (sides == 3) - DrawTriangle(x, y, z, minv, maxv); + DrawTriangle(pts, col, minv, maxv); else - DrawQuad(x, y, z, minv, maxv); + DrawQuad(pts, col, minv, maxv); } else if (split_quads == 1) { @@ -855,46 +942,48 @@ void VisualizationSceneSolution :: PrepareFlat2 () { for (j = 0; j < 3; j++) { - x[j] = pointmat(0, ind[vt[it][j]]); - y[j] = pointmat(1, ind[vt[it][j]]); - z[j] = values(ind[vt[it][j]]); + pts[j][0] = pointmat(0, ind[vt[it][j]]); + pts[j][1] = pointmat(1, ind[vt[it][j]]); + pts[j][2] = col[j] = values(ind[vt[it][j]]); } - DrawTriangle(x, y, z, minv, maxv); + DrawTriangle(pts, col, minv, maxv); } } else { // draw 4 triangles for each quad // (split with both diagonals) - x[2] = y[2] = z[2] = 0.0; + pts[2][0] = pts[2][1] = pts[2][2] = 0.0; for (j = 0; j < 4; j++) { - x[2] += pointmat(0, ind[j]); - y[2] += pointmat(1, ind[j]); - z[2] += values(ind[j]); + pts[2][0] += pointmat(0, ind[j]); + pts[2][1] += pointmat(1, ind[j]); + pts[2][2] += values(ind[j]); } - x[2] *= 0.25; - y[2] *= 0.25; - z[2] *= 0.25; + pts[2][0] *= 0.25; + pts[2][1] *= 0.25; + pts[2][2] *= 0.25; + col[2] = pts[2][2]; for (j = 0; j < 4; j++) { - x[0] = pointmat(0, ind[j]); - y[0] = pointmat(1, ind[j]); - z[0] = values(ind[j]); + pts[0][0] = pointmat(0, ind[j]); + pts[0][1] = pointmat(1, ind[j]); + pts[0][2] = col[0] = values(ind[j]); int l = (j+1)%4; - x[1] = pointmat(0, ind[l]); - y[1] = pointmat(1, ind[l]); - z[1] = values(ind[l]); - DrawTriangle(x, y, z, minv, maxv); + pts[1][0] = pointmat(0, ind[l]); + pts[1][1] = pointmat(1, ind[l]); + pts[1][2] = col[1] = values(ind[l]); + DrawTriangle(pts, col, minv, maxv); } } } +#endif } - glEndList (); + glEndList(); } -void VisualizationSceneSolution :: Prepare () +void VisualizationSceneSolution::Prepare() { switch (shading) { @@ -912,23 +1001,25 @@ void VisualizationSceneSolution :: Prepare () glNewList (displlist, GL_COMPILE); - Set_Material (); + Set_Material(); int ne = mesh -> GetNE(); int nv = mesh -> GetNV(); DenseMatrix pointmat; Array vertices; + double p[4][3], nor[3]; Vector nx(nv); Vector ny(nv); Vector nz(nv); - for (int d = 0; d < mesh -> attributes.Size(); d++){ + for (int d = 0; d < mesh -> attributes.Size(); d++) + { if (!el_attr_to_show[mesh -> attributes[d]-1]) continue; - nx=0; - ny=0; - nz=0; + nx = 0.; + ny = 0.; + nz = 0.; for (i = 0; i < ne; i++) if (mesh -> GetAttribute(i) == mesh -> attributes[d]) @@ -936,24 +1027,25 @@ void VisualizationSceneSolution :: Prepare () mesh->GetPointMatrix (i, pointmat); mesh->GetElementVertices (i, vertices); - double v10[] = { pointmat(0,1)-pointmat(0,0), - pointmat(1,1)-pointmat(1,0), - (*sol)(vertices[1])-(*sol)(vertices[0]) }; - double v21[] = { pointmat(0,2)-pointmat(0,1), - pointmat(1,2)-pointmat(1,1), - (*sol)(vertices[2])-(*sol)(vertices[1]) }; - double norm[] = { v10[1]*v21[2]-v10[2]*v21[1], - v10[2]*v21[0]-v10[0]*v21[2], - v10[0]*v21[1]-v10[1]*v21[0] }; - double rlen = 1.0/sqrt(norm[0]*norm[0]+norm[1]*norm[1]+ - norm[2]*norm[2]); - for (j = 0; j < pointmat.Size(); j++) { - nx(vertices[j]) += norm[0]*rlen; - ny(vertices[j]) += norm[1]*rlen; - nz(vertices[j]) += norm[2]*rlen; + p[j][0] = pointmat(0, j); + p[j][1] = pointmat(1, j); + p[j][2] = (*sol)(vertices[j]); } + + if (pointmat.Width() == 3) + j = Compute3DUnitNormal(p[0], p[1], p[2], nor); + else + j = Compute3DUnitNormal(p[0], p[1], p[2], p[3], nor); + + if (j == 0) + for (j = 0; j < pointmat.Size(); j++) + { + nx(vertices[j]) += nor[0]; + ny(vertices[j]) += nor[1]; + nz(vertices[j]) += nor[2]; + } } for (i = 0; i < ne; i++) @@ -974,17 +1066,15 @@ void VisualizationSceneSolution :: Prepare () for (j = 0; j < pointmat.Size(); j++) { - MySetColor ( (*sol)(vertices[j]) , minv , maxv); - glNormal3d (nx(vertices[j]),ny(vertices[j]),nz(vertices[j])); - glVertex3d (pointmat(0, j), - pointmat(1, j), - (*sol)(vertices[j])); + MySetColor((*sol)(vertices[j]), minv, maxv); + glNormal3d(nx(vertices[j]), ny(vertices[j]), nz(vertices[j])); + glVertex3d(pointmat(0, j), pointmat(1, j), (*sol)(vertices[j])); } - glEnd (); + glEnd(); } } - glEndList (); + glEndList(); } void VisualizationSceneSolution::PrepareLevelCurves() @@ -1015,7 +1105,7 @@ void VisualizationSceneSolution::PrepareLevelCurves() DrawPolygonLevelLines (point[0], vertices.Size(), level); } - glEndList (); + glEndList(); } void VisualizationSceneSolution::DrawLevelCurves( @@ -1115,14 +1205,15 @@ void VisualizationSceneSolution::PrepareLevelCurves2() DrawLevelCurves(RG, pointmat, values, sides, level); } - glEndList (); + glEndList(); } -void VisualizationSceneSolution::PrepareLines(){ +void VisualizationSceneSolution::PrepareLines() +{ if (shading == 2) { - // PrepareLines2 (); - PrepareLines3 (); + // PrepareLines2(); + PrepareLines3(); return; } @@ -1132,7 +1223,7 @@ void VisualizationSceneSolution::PrepareLines(){ glNewList(linelist, GL_COMPILE); - // Set_Black_Material (); + // Set_Black_Material(); // glColor3f (0, 0, 0); glPolygonMode (GL_FRONT_AND_BACK, GL_LINE); @@ -1155,10 +1246,10 @@ void VisualizationSceneSolution::PrepareLines(){ for (j = 0; j < pointmat.Size(); j++) glVertex3d (pointmat(0, j), pointmat(1, j), (*sol)(vertices[j])); - glEnd (); + glEnd(); } - glEndList (); + glEndList(); } void VisualizationSceneSolution::PrepareLines2() @@ -1170,7 +1261,7 @@ void VisualizationSceneSolution::PrepareLines2() glNewList(linelist, GL_COMPILE); - // Set_Black_Material (); + // Set_Black_Material(); // glColor3f (0, 0, 0); glPolygonMode (GL_FRONT_AND_BACK, GL_LINE); @@ -1201,11 +1292,11 @@ void VisualizationSceneSolution::PrepareLines2() glVertex3d (pointmat(0, RG[sides*k+j]), pointmat(1, RG[sides*k+j]), values(RG[sides*k+j])); - glEnd (); + glEnd(); } } - glEndList (); + glEndList(); } void VisualizationSceneSolution::PrepareLines3() @@ -1217,7 +1308,7 @@ void VisualizationSceneSolution::PrepareLines3() glNewList(linelist, GL_COMPILE); - // Set_Black_Material (); + // Set_Black_Material(); // glColor3f (0, 0, 0); // glPolygonMode (GL_FRONT_AND_BACK, GL_LINE); @@ -1239,10 +1330,10 @@ void VisualizationSceneSolution::PrepareLines3() pointmat(1, RE[2*k+1]), values(RE[2*k+1])); } - glEnd (); + glEnd(); } - glEndList (); + glEndList(); } void VisualizationSceneSolution::UpdateValueRange() @@ -1269,10 +1360,10 @@ void VisualizationSceneSolution::PrepareBoundary() glBegin (GL_LINE_LOOP); for (j = 0; j < pointmat.Size(); j++) glVertex3d (pointmat(0, j), pointmat(1, j), (*sol)(vertices[j])); - glEnd (); + glEnd(); } - glEndList (); + glEndList(); } void VisualizationSceneSolution::PrepareCP() @@ -1394,10 +1485,11 @@ void VisualizationSceneSolution::DrawCPLine( } } -void VisualizationSceneSolution :: Draw (){ +void VisualizationSceneSolution::Draw() +{ glEnable(GL_DEPTH_TEST); - Set_Background (); + Set_Background(); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); // model transformation @@ -1422,10 +1514,11 @@ void VisualizationSceneSolution :: Draw (){ else DrawColorBar(minv,maxv); - Set_Black_Material (); + Set_Black_Material(); // draw axes - if (drawaxes){ + if (drawaxes) + { glCallList(axeslist); DrawCoordinateCross(); } @@ -1444,10 +1537,11 @@ void VisualizationSceneSolution :: Draw (){ glCallList(bdrlist); // draw lines - if ( drawmesh == 1 ) + if (drawmesh == 1) glCallList(linelist); else - if (drawmesh == 2 ){ + if (drawmesh == 2) + { // glLineWidth(1.0); glCallList(lcurvelist); } @@ -1455,7 +1549,7 @@ void VisualizationSceneSolution :: Draw (){ glEnable(GL_LIGHTING); if (MatAlpha < 1.0) - Set_Transparency (); + Set_Transparency(); // draw elements glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); @@ -1472,7 +1566,7 @@ void VisualizationSceneSolution :: Draw (){ } if (MatAlpha < 1.0) - Remove_Transparency (); + Remove_Transparency(); glFlush(); glXSwapBuffers (auxXDisplay(), auxXWindow()); diff --git a/lib/vssolution.hpp b/lib/vssolution.hpp index 8ab9aee..405620f 100644 --- a/lib/vssolution.hpp +++ b/lib/vssolution.hpp @@ -34,12 +34,12 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData void GetRefinedDetJ(int i, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr); - /// Shrink the set of points towards their center of gravity - void ShrinkPoints(DenseMatrix &pointmat); - // redefined for vector solution virtual void GetRefinedValues(int i, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr); + virtual int GetRefinedValuesAndNormals(int i, const IntegrationRule &ir, + Vector &vals, DenseMatrix &tr, + DenseMatrix &normals); void DrawLevelCurves(Array &RG, DenseMatrix &pointmat, Vector &values, int sides, Array &lvl, @@ -47,7 +47,6 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData public: int shading, TimesToRefine, EdgeRefineFactor; - double shrink; int attr_to_show; Array el_attr_to_show; @@ -86,10 +85,10 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData virtual void Draw(); - void ToggleDrawBdr () + void ToggleDrawBdr() { drawbdr = !drawbdr; } - void ToggleDrawElems () + void ToggleDrawElems() { drawelems = (drawelems+3)%4; if (drawelems != 0 && shading == 2) @@ -104,8 +103,7 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData } } - void ToggleDrawMesh () - { drawmesh = (drawmesh+1)%3; } + void ToggleDrawMesh() { drawmesh = (drawmesh+1)%3; } virtual void SetShading(int); void ToggleShading(); @@ -116,4 +114,14 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData }; +void DrawTriangle(const double pts[][3], const double cv[], + const double minv, const double maxv); + +void DrawQuad(const double pts[][3], const double cv[], + const double minv, const double maxv); + +void DrawPatch(const DenseMatrix &pts, Vector &vals, DenseMatrix &normals, + const int n, const Array &ind, const double minv, + const double maxv, const int normals_opt = 0); + #endif diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 9d74044..ec369ae 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -18,16 +18,13 @@ #include "mfem.hpp" #include "visual.hpp" -void Set_Material(); -void Set_Black_Material(); -VisualizationSceneSolution3d * vssol3d; -extern VisualizationScene * locscene; -extern VisualizationSceneScalarData * vsdata; +VisualizationSceneSolution3d *vssol3d; // Definitions of some more keys -static void Solution3dKeyHPressed (){ +static void Solution3dKeyHPressed() +{ cout << endl << "+------------------------------------+" << endl << "| Keys |" << endl @@ -69,13 +66,13 @@ static void Solution3dKeyHPressed (){ << "+------------------------------------+" << endl << "| F1 - X window info and keystrokes |" << endl << "| F2 - Update colors, etc. |" << endl - << "| F3/F4 - Shrink/Zoom subdomains |" << endl + << "| F3/F4 - Shrink/Zoom bdr elements |" << endl << "| F5 - Set level lines |" << endl << "| F6 - Palete options |" << endl << "| F7 - Manually set min/max value |" << endl << "| F8 - List of subdomains to show |" << endl << "| F9/F10 - Walk through subdomains |" << endl - << "| F11/F12 - Shrink/Zoom elements |" << endl + << "| F11/F12 - Shrink/Zoom subdomains |" << endl << "+------------------------------------+" << endl << "| Keypad |" << endl << "+------------------------------------+" << endl @@ -100,7 +97,7 @@ static void Solution3dKeyHPressed (){ << "+------------------------------------+" << endl; } -static void KeyIPressed () +static void KeyIPressed() { vssol3d -> ToggleCuttingPlane(); SendExposeEvent(); @@ -110,19 +107,23 @@ void VisualizationSceneSolution3d::CPPrepare() { PrepareCuttingPlane(); PrepareCuttingPlaneLines(); +} + +void VisualizationSceneSolution3d::CPMoved() +{ + CPPrepare(); if (cplane == 2) { Prepare(); PrepareLines(); } - SendExposeEvent(); } static void KeyxPressed() { vssol3d -> CuttingPlane -> IncreasePhi(); vssol3d -> FindNodePos(); - vssol3d -> CPPrepare(); + vssol3d -> CPMoved(); SendExposeEvent(); } @@ -130,7 +131,7 @@ static void KeyXPressed() { vssol3d -> CuttingPlane -> DecreasePhi(); vssol3d -> FindNodePos(); - vssol3d -> CPPrepare(); + vssol3d -> CPMoved(); SendExposeEvent(); } @@ -138,7 +139,7 @@ static void KeyyPressed() { vssol3d -> CuttingPlane -> IncreaseTheta(); vssol3d -> FindNodePos(); - vssol3d -> CPPrepare(); + vssol3d -> CPMoved(); SendExposeEvent(); } @@ -146,7 +147,7 @@ static void KeyYPressed() { vssol3d -> CuttingPlane -> DecreaseTheta(); vssol3d -> FindNodePos(); - vssol3d -> CPPrepare(); + vssol3d -> CPMoved(); SendExposeEvent(); } @@ -154,7 +155,7 @@ static void KeyzPressed() { vssol3d -> CuttingPlane -> IncreaseDistance(); vssol3d -> FindNodePos(); - vssol3d -> CPPrepare(); + vssol3d -> CPMoved(); SendExposeEvent(); } @@ -162,7 +163,7 @@ static void KeyZPressed() { vssol3d -> CuttingPlane -> DecreaseDistance(); vssol3d -> FindNodePos(); - vssol3d -> CPPrepare(); + vssol3d -> CPMoved(); SendExposeEvent(); } @@ -204,7 +205,7 @@ static void KeyoPressed() << endl; if (vssol3d -> GetShading() == 2) { - locscene -> Prepare(); + vssol3d -> Prepare(); vssol3d -> PrepareLines(); vssol3d -> CPPrepare(); SendExposeEvent(); @@ -220,7 +221,7 @@ static void KeyOPressed() << endl; if (vssol3d -> GetShading() == 2) { - locscene -> Prepare(); + vssol3d -> Prepare(); vssol3d -> PrepareLines(); vssol3d -> CPPrepare(); SendExposeEvent(); @@ -236,7 +237,7 @@ static void KeywPressed() vssol3d -> FaceShiftScale = 0.0; cout << "New Shift Scale: " << vssol3d -> FaceShiftScale << endl; - locscene -> Prepare(); + vssol3d -> Prepare(); vssol3d -> PrepareLines(); vssol3d -> CPPrepare(); SendExposeEvent(); @@ -251,7 +252,7 @@ static void KeyWPressed() vssol3d -> FaceShiftScale = 0.0; cout << "New Shift Scale: " << vssol3d -> FaceShiftScale << endl; - locscene -> Prepare(); + vssol3d -> Prepare(); vssol3d -> PrepareLines(); vssol3d -> CPPrepare(); SendExposeEvent(); @@ -290,65 +291,91 @@ void ToggleMagicKey() void KeyF3Pressed() { - (vsdata -> GetMesh()) -> ScaleSubdomains(0.9); - if (magic_key_pressed) - locscene -> Scale(1.11111111111111111111111); - vssol3d -> PrepareLines(); - vssol3d -> Prepare(); - SendExposeEvent(); + if (vssol3d->GetShading() == 2) + { + if (vssol3d->bdrc.Width() == 0) + vssol3d->ComputeBdrAttrCenter(); + vssol3d->shrink *= 0.9; + if (magic_key_pressed) + vssol3d -> Scale(1.11111111111111111111111); + vssol3d->Prepare(); + vssol3d->PrepareLines(); + SendExposeEvent(); + } } void KeyF4Pressed() { - (vsdata -> GetMesh()) -> ScaleSubdomains(1.11111111111111111111111); - if (magic_key_pressed) - locscene -> Scale(0.9); - vssol3d -> PrepareLines(); - vssol3d -> Prepare(); - SendExposeEvent(); + if (vssol3d->GetShading() == 2) + { + if (vssol3d->bdrc.Width() == 0) + vssol3d->ComputeBdrAttrCenter(); + vssol3d->shrink *= 1.11111111111111111111111; + if (magic_key_pressed) + vssol3d -> Scale(0.9); + vssol3d->Prepare(); + vssol3d->PrepareLines(); + SendExposeEvent(); + } } void KeyF11Pressed() { - (vsdata -> GetMesh()) -> ScaleElements(0.9); - vssol3d -> PrepareLines(); - vssol3d -> Prepare(); - SendExposeEvent(); + if (vssol3d->GetShading() == 2) + { + if (vssol3d->matc.Width() == 0) + vssol3d->ComputeElemAttrCenter(); + vssol3d->shrinkmat *= 0.9; + if (magic_key_pressed) + vssol3d -> Scale(1.11111111111111111111111); + vssol3d->Prepare(); + vssol3d->PrepareLines(); + SendExposeEvent(); + } } void KeyF12Pressed() { - (vsdata -> GetMesh()) -> ScaleElements(1.11111111111111111111111); - vssol3d -> PrepareLines(); - vssol3d -> Prepare(); - SendExposeEvent(); + if (vssol3d->GetShading() == 2) + { + if (vssol3d->matc.Width() == 0) + vssol3d->ComputeElemAttrCenter(); + vssol3d->shrinkmat *= 1.11111111111111111111111; + if (magic_key_pressed) + vssol3d -> Scale(0.9); + vssol3d->Prepare(); + vssol3d->PrepareLines(); + SendExposeEvent(); + } } static void KeyF8Pressed() { + const Array &attr_list = vssol3d->GetMesh()->bdr_attributes; + Array &attr_marker = vssol3d->bdr_attr_to_show; int attr; cout << "Bdr attributes ON: "; - for (attr = 0; attr < vssol3d -> bdr_attr_to_show.Size(); attr++) - if (vssol3d -> bdr_attr_to_show[attr]) - cout << " " << attr; + for (int i = 0; i < attr_list.Size(); i++) + if (attr_marker[attr_list[i]-1]) + cout << " " << attr_list[i]; cout << endl; cout << "Bdr attribute to toggle : " << flush; cin >> attr; - if (attr < 0 ) + if (attr < 1) { - for (attr = 0; attr < vssol3d -> bdr_attr_to_show.Size(); attr++) - vssol3d -> bdr_attr_to_show[attr] = 0; + cout << "Hiding all bdr attributes." << endl; + attr_marker = 0; } - else if (attr >= vssol3d -> bdr_attr_to_show.Size()) + else if (attr > attr_marker.Size()) { - for (attr = 0; attr < vssol3d -> bdr_attr_to_show.Size(); attr++) - vssol3d -> bdr_attr_to_show[attr] = 1; + cout << "Showing all bdr attributes." << endl; + attr_marker = 1; } else { - vssol3d -> bdr_attr_to_show[attr] = !vssol3d -> bdr_attr_to_show[attr]; + attr_marker[attr-1] = !attr_marker[attr-1]; } vssol3d -> PrepareLines(); vssol3d -> Prepare(); @@ -357,19 +384,24 @@ static void KeyF8Pressed() static void KeyF9Pressed() { - int attr; + const Array &attr_list = vssol3d->GetMesh()->bdr_attributes; + Array &attr_marker = vssol3d->bdr_attr_to_show; + int attr, n, j; - if (vssol3d -> attr_to_show == -1) { - for (attr = 0; attr < vssol3d -> bdr_attr_to_show.Size(); attr++) - vssol3d -> bdr_attr_to_show[attr] = 0; - vssol3d -> attr_to_show = 0; - } else - vssol3d -> bdr_attr_to_show[vssol3d -> attr_to_show++] = 0; - - if (vssol3d -> attr_to_show >= vssol3d -> bdr_attr_to_show.Size()) - vssol3d -> attr_to_show = 0; - vssol3d -> bdr_attr_to_show[vssol3d -> attr_to_show] = 1; - cout << "Showing bdr attribute " << vssol3d -> attr_to_show << endl; + if (attr_list.Size() == 0) + return; + n = 0; + for (int i = 0; i < attr_list.Size(); i++) + if (attr_marker[attr_list[i]-1]) + j = i, n++; + if (n == 1) + j = (j + 1) % attr_list.Size(); + else + j = 0; + attr = attr_list[j]; + attr_marker = 0; + attr_marker[attr-1] = 1; + cout << "Showing bdr attribute " << attr << endl; vssol3d -> PrepareLines(); vssol3d -> Prepare(); SendExposeEvent(); @@ -377,91 +409,101 @@ static void KeyF9Pressed() static void KeyF10Pressed() { - int attr; + const Array &attr_list = vssol3d->GetMesh()->bdr_attributes; + Array &attr_marker = vssol3d->bdr_attr_to_show; + int attr, n, j; - if (vssol3d -> attr_to_show == -1) { - for (attr = 0; attr < vssol3d -> bdr_attr_to_show.Size(); attr++) - vssol3d -> bdr_attr_to_show[attr] = 0; - } else - vssol3d -> bdr_attr_to_show[vssol3d -> attr_to_show--] = 0; - - if (vssol3d -> attr_to_show < 0) - vssol3d -> attr_to_show = vssol3d -> bdr_attr_to_show.Size()-1; - vssol3d -> bdr_attr_to_show[vssol3d -> attr_to_show] = 1; - cout << "Showing bdr attribute " << vssol3d -> attr_to_show << endl; + if (attr_list.Size() == 0) + return; + n = 0; + for (int i = 0; i < attr_list.Size(); i++) + if (attr_marker[attr_list[i]-1]) + j = i, n++; + if (n == 1) + j = (j + attr_list.Size() - 1) % attr_list.Size(); + else + j = attr_list.Size() - 1; + attr = attr_list[j]; + attr_marker = 0; + attr_marker[attr-1] = 1; + cout << "Showing bdr attribute " << attr << endl; vssol3d -> PrepareLines(); vssol3d -> Prepare(); SendExposeEvent(); } -double Distance (double p1[], double p2[]) +// 'v' must define three vectors that add up to the zero vector +// that is v[0], v[1], and v[2] are the sides of a triangle +int UnitCrossProd(double v[][3], double nor[]) { - double v[] = { p2[0]-p1[0], p2[1]-p1[1], p2[2]-p1[2] }; + // normalize the three vectors + for (int i = 0; i < 3; i++) + if (Normalize(v[i])) + return 1; + + // find the pair that forms an angle closest to pi/2, i.e. having + // the longest cross product: + double cp[3][3], max_a = 0.; + int k = 0; + for (int i = 0; i < 3; i++) + { + CrossProd(v[(i+1)%3], v[(i+2)%3], cp[i]); + double a = sqrt(InnerProd(cp[i], cp[i])); + if (max_a < a) + max_a = a, k = i; + } + if (max_a == 0.) + return 1; + for (int i = 0; i < 3; i++) + nor[i] = cp[k][i] / max_a; - return sqrt (v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); + return 0; } -double InnerProduct (double p1[], double p2[], double p3[]) +int Compute3DUnitNormal(const double p1[], const double p2[], + const double p3[], double nor[]) { - double v1[] = { p2[0]-p1[0], p2[1]-p1[1], p2[2]-p1[2] }; - double v2[] = { p3[0]-p1[0], p3[1]-p1[1], p3[2]-p1[2] }; + double v[3][3]; - return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]; + for (int i = 0; i < 3; i++) + { + v[0][i] = p2[i] - p1[i]; + v[1][i] = p3[i] - p2[i]; + v[2][i] = p1[i] - p3[i]; + } + + return UnitCrossProd(v, nor); } -int Compute3DUnitNormal (double p1[], double p2[], double p3[], double nor[]) +int Compute3DUnitNormal (const double p1[], const double p2[], + const double p3[], const double p4[], double nor[]) { - int err = 0; - double d1 = Distance (p2, p3); - double d2 = Distance (p1, p3); - double d3 = Distance (p1, p2); - double *pp[3], pr, dmin; - if (d1 < d2) - if (d1 < d3) - pp[0] = p1, pp[1] = p2, pp[2] = p3, pr = d2*d3, dmin = d1; - else - pp[0] = p3, pp[1] = p1, pp[2] = p2, pr = d1*d2, dmin = d3; - else - if (d2 < d3) - pp[0] = p2, pp[1] = p3, pp[2] = p1, pr = d1*d3, dmin = d2; - else - pp[0] = p3, pp[1] = p1, pp[2] = p2, pr = d1*d2, dmin = d3; - double inpr = InnerProduct (pp[0], pp[1], pp[2]); - if (dmin == 0.0 || inpr >= (1.-1.e-9) * pr) - err += 2; - double v1[] = { pp[1][0]-pp[0][0], pp[1][1]-pp[0][1], pp[1][2]-pp[0][2] }; - double v2[] = { pp[2][0]-pp[0][0], pp[2][1]-pp[0][1], pp[2][2]-pp[0][2] }; - double n[] = { v1[1] * v2[2] - v1[2] * v2[1], - v1[2] * v2[0] - v1[0] * v2[2], - v1[0] * v2[1] - v1[1] * v2[0] }; - double rlen = 1.0 / sqrt (n[0]*n[0]+n[1]*n[1]+n[2]*n[2]); - - if (!finite(rlen)) - rlen = 1.0, err += 1; - - nor[0] = rlen * n[0]; - nor[1] = rlen * n[1]; - nor[2] = rlen * n[2]; + double v[3][3]; - return err; -} + for (int i = 0; i < 3; i++) + { + // cross product of the two diagonals: + /* + v[0][i] = p3[i] - p1[i]; + v[1][i] = p4[i] - p2[i]; + v[2][i] = (p1[i] + p2[i]) - (p3[i] + p4[i]); + */ -int Compute3DUnitNormal (double p1[], double p2[], double p3[], double p4[], - double nor[]) -{ - if (Compute3DUnitNormal(p1, p2, p3, nor) == 0) - return 0; - if (Compute3DUnitNormal(p1, p2, p4, nor) == 0) - return 0; - if (Compute3DUnitNormal(p1, p3, p4, nor) == 0) - return 0; - return Compute3DUnitNormal(p2, p3, p4, nor); + // cross product of the two vectors connecting the midpoints of the + // two pairs of opposing sides; this gives the normal vector in the + // midpoint of the quad: + v[0][i] = 0.5 * ((p2[i] + p3[i]) - (p1[i] + p4[i])); + v[1][i] = 0.5 * ((p4[i] + p3[i]) - (p1[i] + p2[i])); + v[2][i] = p1[i] - p3[i]; + } + + return UnitCrossProd(v, nor); } -VisualizationSceneSolution3d::VisualizationSceneSolution3d(){} +VisualizationSceneSolution3d::VisualizationSceneSolution3d() +{} -VisualizationSceneSolution3d -::VisualizationSceneSolution3d ( Mesh & m, Vector & s ) +VisualizationSceneSolution3d::VisualizationSceneSolution3d(Mesh &m, Vector &s) { mesh = &m; sol = &s; @@ -487,23 +529,24 @@ void VisualizationSceneSolution3d::Init() drawmesh = 0; scaling = 0; + shrink = 1.0; + shrinkmat = 1.0; + bdrc.SetSize(3,0); + matc.SetSize(3,0); + TimesToRefine = 1; FaceShiftScale = 0.0; - attr_to_show = -1; - minv = sol->Min(); maxv = sol->Max(); - int i; bdr_attr_to_show.SetSize (mesh->bdr_attributes.Max()); - for (i = 0; i < bdr_attr_to_show.Size(); i++) - bdr_attr_to_show[i] = 1; + bdr_attr_to_show = 1; - VisualizationSceneScalarData :: Init(); + VisualizationSceneScalarData::Init(); - scaling = 0; // - SetNewScalingFromBox (); // No scaling for 3D + scaling = 0; + SetNewScalingFromBox(); // No scaling for 3D node_pos = new double[mesh->GetNV()]; @@ -518,7 +561,6 @@ void VisualizationSceneSolution3d::Init() // This implements a) SetLevelLines (minv, maxv, 15); // 15 level curves - nlevels = 1; FindNodePos(); @@ -584,7 +626,8 @@ void VisualizationSceneSolution3d::Init() PrepareLevelSurf(); } -VisualizationSceneSolution3d::~VisualizationSceneSolution3d (){ +VisualizationSceneSolution3d::~VisualizationSceneSolution3d() +{ glDeleteLists (displlist, 1); glDeleteLists (linelist, 1); glDeleteLists (cplanelist, 1); @@ -593,6 +636,38 @@ VisualizationSceneSolution3d::~VisualizationSceneSolution3d (){ delete [] node_pos; } +void VisualizationSceneSolution3d::NewMeshAndSolution( + Mesh *new_m, Vector *new_sol, GridFunction *new_u, int rescale) +{ + if (mesh->GetNV() != new_m->GetNV()) + { + delete [] node_pos; + node_pos = new double[new_m->GetNV()]; + } + mesh = new_m; + sol = new_sol; + GridF = new_u; + FindNodePos(); + + if (rescale == 1) + { + FindNewBox(); + PrepareAxes(); + minv = sol->Min(); + maxv = sol->Max(); + } + else if (rescale == 2) + { + minv = sol->Min(); + maxv = sol->Max(); + } + + Prepare(); + PrepareLines(); + CPPrepare(); + PrepareLevelSurf(); +} + void VisualizationSceneSolution3d::SetShading(int s) { if (shading == s || s < 0) @@ -608,6 +683,10 @@ void VisualizationSceneSolution3d::SetShading(int s) PrepareLines(); CPPrepare(); } + + static const char *shading_type[3] = + {"flat", "smooth", "non-conforming (with subdivision)"}; + cout << "Shading type : " << shading_type[shading] << endl; } void VisualizationSceneSolution3d::ToggleShading() @@ -656,14 +735,14 @@ void VisualizationSceneSolution3d::FindNewBox() if (coord[2] > z[1]) z[1] = coord[2]; } - SetNewScalingFromBox (); + SetNewScalingFromBox(); } -void VisualizationSceneSolution3d :: FindNodePos() +void VisualizationSceneSolution3d::FindNodePos() { int i, nnodes = mesh -> GetNV(); - for(i = 0; i < nnodes; i++) + for (i = 0; i < nnodes; i++) node_pos[i] = CuttingPlane -> Transform (mesh -> GetVertex (i)); } @@ -681,22 +760,22 @@ void VisualizationSceneSolution3d::ToggleCuttingPlane() cplane = (cplane+1)%3; if (cplane) - CPPrepare (); - else + CPPrepare(); + if (cplane == 0 || cplane == 2) { Prepare(); PrepareLines(); } } -void VisualizationSceneSolution3d::ToggleCPDrawElems () +void VisualizationSceneSolution3d::ToggleCPDrawElems() { cp_drawelems = 1-cp_drawelems; if (cp_drawelems) PrepareCuttingPlane(); } -void VisualizationSceneSolution3d::ToggleCPDrawMesh () +void VisualizationSceneSolution3d::ToggleCPDrawMesh() { if (cplane == 1) cp_drawmesh = (cp_drawmesh+1)%3; @@ -706,25 +785,76 @@ void VisualizationSceneSolution3d::ToggleCPDrawMesh () PrepareCuttingPlaneLines(); } -void VisualizationSceneSolution3d::MoveLevelSurf (int move) +void VisualizationSceneSolution3d::MoveLevelSurf(int move) { drawlsurf += move; if (drawlsurf < 0) drawlsurf = 0; if (drawlsurf > 49) drawlsurf = 49; - PrepareLevelSurf (); + PrepareLevelSurf(); } -void VisualizationSceneSolution3d::NumberOfLevelSurf (int c) +void VisualizationSceneSolution3d::NumberOfLevelSurf(int c) { nlevels += c; if (nlevels < 1) nlevels = 1; - PrepareLevelSurf (); + PrepareLevelSurf(); +} + +int Normalize(DenseMatrix &normals) +{ + int err = 0; + for (int i = 0; i < normals.Width(); i++) + err += Normalize(&normals(0, i)); + return err; } -void VisualizationSceneSolution3d::DrawRefinedSurf ( +void VisualizationSceneSolution3d::GetFaceNormals( + const int FaceNo, const int side, const IntegrationRule &ir, + DenseMatrix &normals) +{ + // match the way GridFunction::GetFaceValues works + double aJInv[9], alnor[3]; + DenseMatrix JInv(aJInv, 3, 3); + Vector lnor(alnor, 3), nr; + IntegrationRule eir(ir.GetNPoints()); + FaceElementTransformations *Tr = + mesh->GetFaceElementTransformations(FaceNo); + normals.SetSize(3, ir.GetNPoints()); + ElementTransformation *ETr; + IntegrationPointTransformation *LTr; + if (side == 0) + { + ETr = Tr->Elem1; + LTr = &Tr->Loc1; + } + else + { + ETr = Tr->Elem2; + LTr = &Tr->Loc2; + } + LTr->Transform(ir, eir); + for (int i = 0; i < normals.Width(); i++) + { + LTr->Transf.SetIntPoint(&ir.IntPoint(i)); + const DenseMatrix &LJac = LTr->Transf.Jacobian(); + lnor(0) = LJac(1,0) * LJac(2,1) - LJac(2,0) * LJac(1,1); + lnor(1) = LJac(2,0) * LJac(0,1) - LJac(0,0) * LJac(2,1); + lnor(2) = LJac(0,0) * LJac(1,1) - LJac(1,0) * LJac(0,1); + ETr->SetIntPoint(&eir.IntPoint(i)); + const DenseMatrix &Jac = ETr->Jacobian(); + CalcInverse(Jac, JInv); + normals.GetColumnReference(i, nr); + JInv.MultTranspose(lnor, nr); + } + if (side) + normals *= -1.; + JInv.ClearExternalData(); +} + +void VisualizationSceneSolution3d::DrawRefinedSurf( int n, double *points, int elem, int func, int part) { int i, j; @@ -785,7 +915,7 @@ void VisualizationSceneSolution3d::DrawRefinedSurf ( } } -void VisualizationSceneSolution3d::LiftRefinedSurf ( +void VisualizationSceneSolution3d::LiftRefinedSurf( int n, DenseMatrix &pointmat, Vector &values, int *RG) { int i, j; @@ -834,7 +964,7 @@ void VisualizationSceneSolution3d::LiftRefinedSurf ( } } -void VisualizationSceneSolution3d::DrawRefinedSurf ( +void VisualizationSceneSolution3d::DrawRefinedSurf( int n, DenseMatrix &pointmat, Vector &values, Array &RefGeoms) { double norm[3], pts[4][3]; @@ -859,7 +989,7 @@ void VisualizationSceneSolution3d::DrawRefinedSurf ( MySetColor (values(RG[j]), minv, maxv); glVertex3dv (pts[j]); } - glEnd (); + glEnd(); } /* else @@ -869,7 +999,7 @@ void VisualizationSceneSolution3d::DrawRefinedSurf ( } } -void VisualizationSceneSolution3d::DrawRefinedSurfEdges ( +void VisualizationSceneSolution3d::DrawRefinedSurfEdges( int n, DenseMatrix &pointmat, Vector &values, Array &RefEdges, int part) { @@ -892,10 +1022,10 @@ void VisualizationSceneSolution3d::DrawRefinedSurfEdges ( pointmat(2, RE)); } if (part != 0) - glEnd (); + glEnd(); } -void VisualizationSceneSolution3d::DrawRefinedSurfLevelLines ( +void VisualizationSceneSolution3d::DrawRefinedSurfLevelLines( int n, DenseMatrix &pointmat, Vector &values, Array &RefGeoms) { int j, k; @@ -916,17 +1046,19 @@ void VisualizationSceneSolution3d::DrawRefinedSurfLevelLines ( } } -void VisualizationSceneSolution3d :: PrepareFlat (){ +void VisualizationSceneSolution3d::PrepareFlat() +{ int i, j; glNewList (displlist, GL_COMPILE); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); - Set_Material (); + Set_Material(); int ne = mesh -> GetNBE(); DenseMatrix pointmat; Array vertices; + double p[4][3], c[4]; for (i = 0; i < ne; i++) { @@ -944,58 +1076,39 @@ void VisualizationSceneSolution3d :: PrepareFlat (){ continue; // with the next boundary element } - switch (mesh->GetBdrElementType(i)) - { - case Element::TRIANGLE: - glBegin (GL_TRIANGLES); - break; - - case Element::QUADRILATERAL: - glBegin (GL_QUADS); - break; - } mesh->GetBdrPointMatrix (i, pointmat); - double v10[] = { pointmat(0,1)-pointmat(0,0), - pointmat(1,1)-pointmat(1,0), - pointmat(2,1)-pointmat(2,0) }; - double v21[] = { pointmat(0,2)-pointmat(0,1), - pointmat(1,2)-pointmat(1,1), - pointmat(2,2)-pointmat(2,1) }; - - double norm[] = { v10[1]*v21[2]-v10[2]*v21[1], - v10[2]*v21[0]-v10[0]*v21[2], - v10[0]*v21[1]-v10[1]*v21[0] }; - double rlen = 1.0/sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]); - - glNormal3d (norm[0]*rlen, norm[1]*rlen, norm[2]*rlen); - - for (j = 0; j < pointmat.Size(); j++) { - MySetColor ( (*sol)(vertices[j]) , minv, maxv ); - glVertex3d (pointmat(0, j), - pointmat(1, j), - pointmat(2, j)); + for (j = 0; j < pointmat.Width(); j++) + { + p[j][0] = pointmat(0, j); + p[j][1] = pointmat(1, j); + p[j][2] = pointmat(2, j); + c[j] = (*sol)(vertices[j]); } - glEnd (); + if (j == 3) + DrawTriangle(p, c, minv, maxv); + else + DrawQuad(p, c, minv, maxv); } - glEndList (); + glEndList(); } void VisualizationSceneSolution3d::PrepareFlat2() { - int i, j, k, fn, fo, ft, di; + int i, j, k, fn, fo, ft, di, have_normals; double bbox_diam, vmin, vmax, mm; glNewList (displlist, GL_COMPILE); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); - Set_Material (); + Set_Material(); int nbe = mesh -> GetNBE(); - DenseMatrix pointmat; + DenseMatrix pointmat, normals; Vector values; RefinedGeometry * RefG; Array vertices; + double norm[3]; bbox_diam = sqrt ( (x[1]-x[0])*(x[1]-x[0]) + (y[1]-y[0])*(y[1]-y[0]) + @@ -1027,6 +1140,9 @@ void VisualizationSceneSolution3d::PrepareFlat2() // this assumes the interior boundary faces are properly oriented ... di = fo % 2; GridF -> GetFaceValues (fn, di, RefG->RefPts, values, pointmat); + GetFaceNormals(fn, di, RefG->RefPts, normals); + have_normals = 1; + ShrinkPoints3D(pointmat, i, fn, fo); if (ft) { vmin = values.Min(); @@ -1052,89 +1168,43 @@ void VisualizationSceneSolution3d::PrepareFlat2() sides = 4; break; } - int *RG = &(RefG->RefGeoms[0]); - double pts[][3] = { { pointmat(0,RG[0]), pointmat(1,RG[0]), - pointmat(2,RG[0]) }, - { pointmat(0,RG[1]), pointmat(1,RG[1]), - pointmat(2,RG[1]) }, - { pointmat(0,RG[2]), pointmat(1,RG[2]), - pointmat(2,RG[2]) } }; - double norm[3]; - if (Compute3DUnitNormal (pts[0], pts[1], pts[2], norm)) - cerr << "WARNING: VisualizationSceneSolution3d::PrepareFlat2()" - << endl; - if (di != 0 && sc != 0.0) - { - norm[0] = -norm[0]; - norm[1] = -norm[1]; - norm[2] = -norm[2]; - } - for (k = 0; k < RefG->RefGeoms.Size()/sides; k++) + // compute an average normal direction for the current face + if (sc != 0.0) { - RG = &(RefG->RefGeoms[k*sides]); - switch (sides) - { - case 3: - glBegin (GL_TRIANGLES); - break; - - case 4: - glBegin (GL_QUADS); - break; - } - - if (sc == 0.0) + for (int i = 0; i < 3; i++) + norm[i] = 0.0; + Normalize(normals); + for (k = 0; k < normals.Width(); k++) + for (int j = 0; j < 3; j++) + norm[j] += normals(j, k); + Normalize(norm); + for (int i = 0; i < pointmat.Width(); i++) { - glNormal3dv (norm); - for (j = 0; j < sides; j++) - { - MySetColor ( values(RG[j]) , minv , maxv ); - glVertex3d (pointmat(0, RG[j]), pointmat(1, RG[j]), - pointmat(2, RG[j])); - } + double val = sc * (values(i) - minv) / (maxv - minv); + for (int j = 0; j < 3; j++) + pointmat(j, i) += val * norm[j]; } - else - { - double nnorm[3]; - for (j = 0; j < 3; j++) - { - double val = sc * (values(RG[j]) - minv) / (maxv - minv); - for (int l = 0; l < 3; l++) - pts[j][l] = pointmat(l, RG[j]) + val * norm[l]; - } - if (Compute3DUnitNormal (pts[0], pts[1], pts[2], nnorm)) - cerr << "WARNING: VisualizationSceneSolution3d::PrepareFlat2()" - << endl; - glNormal3dv (nnorm); - for (j = 0; j < 3; j++) - { - MySetColor ( values(RG[j]) , minv , maxv ); - glVertex3dv (pts[j]); - } - for ( ; j < sides; j++) - { - double val = (values(RG[j]) - minv) / (maxv - minv); - MySetColor ( values(RG[j]) , minv , maxv ); - val *= sc; - glVertex3d (pointmat(0, RG[j])+val*norm[0], - pointmat(1, RG[j])+val*norm[1], - pointmat(2, RG[j])+val*norm[2]); - } - } - glEnd (); + have_normals = 0; } + + have_normals = have_normals ? 2 : 0; + if (di) + have_normals = -1 - have_normals; + DrawPatch(pointmat, values, normals, sides, RefG->RefGeoms, + minv, maxv, have_normals); } - glEndList (); + glEndList(); cout << "VisualizationSceneSolution3d::PrepareFlat2() : [min,max] = [" << vmin << "," << vmax << "]" << endl; } -void VisualizationSceneSolution3d :: Prepare () +void VisualizationSceneSolution3d::Prepare() { int i,j; - switch (shading){ + switch (shading) + { case 0: PrepareFlat(); return; @@ -1147,88 +1217,102 @@ void VisualizationSceneSolution3d :: Prepare () glNewList (displlist, GL_COMPILE); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); - Set_Material (); + Set_Material(); int ne = mesh -> GetNBE(); int nv = mesh -> GetNV(); DenseMatrix pointmat; Array vertices; + double nor[3]; Vector nx(nv); Vector ny(nv); Vector nz(nv); + Table ba_to_be; // boundary_attribute--to--boundary_element + { + Table be_to_ba; + be_to_ba.MakeI(ne); + for (int i = 0; i < ne; i++) + be_to_ba.AddAColumnInRow(i); + be_to_ba.MakeJ(); + for (int i = 0; i < ne; i++) + be_to_ba.AddConnection(i, mesh->GetBdrAttribute(i)-1); + be_to_ba.ShiftUpI(); + + Transpose(be_to_ba, ba_to_be); + } + for (int d = 0; d < mesh -> bdr_attributes.Size(); d++) { - if (!bdr_attr_to_show[mesh -> bdr_attributes[d]-1]) continue; + const int attr = mesh->bdr_attributes[d]-1; - nx=0.; - ny=0.; - nz=0.; + if (!bdr_attr_to_show[attr]) continue; - for (i = 0; i < ne; i++) - if (mesh -> GetBdrAttribute(i) == mesh -> bdr_attributes[d]) - { - mesh->GetBdrPointMatrix (i, pointmat); - mesh->GetBdrElementVertices (i, vertices); - - double v10[] = { pointmat(0,1)-pointmat(0,0), - pointmat(1,1)-pointmat(1,0), - pointmat(2,1)-pointmat(2,0)}; - double v21[] = { pointmat(0,2)-pointmat(0,1), - pointmat(1,2)-pointmat(1,1), - pointmat(2,2)-pointmat(2,1) }; - double norm[] = { v10[1]*v21[2]-v10[2]*v21[1], - v10[2]*v21[0]-v10[0]*v21[2], - v10[0]*v21[1]-v10[1]*v21[0] }; - double rlen = 1.0/sqrt(norm[0]*norm[0]+norm[1]*norm[1]+ - norm[2]*norm[2]); + const int nelem = ba_to_be.RowSize(attr); + const int *elem = ba_to_be.GetRow(attr); + for (i = 0; i < nelem; i++) + { + mesh->GetBdrElementVertices(elem[i], vertices); + for (j = 0; j < vertices.Size(); j++) + nx(vertices[j]) = ny(vertices[j]) = nz(vertices[j]) = 0.; + } + for (i = 0; i < nelem; i++) + { + mesh->GetBdrPointMatrix(elem[i], pointmat); + mesh->GetBdrElementVertices(elem[i], vertices); + + if (pointmat.Width() == 3) + j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1), + &pointmat(0,2), nor); + else + j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1), + &pointmat(0,2), &pointmat(0,3), nor); + if (j == 0) for (j = 0; j < pointmat.Size(); j++) { - nx(vertices[j]) += norm[0]*rlen; - ny(vertices[j]) += norm[1]*rlen; - nz(vertices[j]) += norm[2]*rlen; + nx(vertices[j]) += nor[0]; + ny(vertices[j]) += nor[1]; + nz(vertices[j]) += nor[2]; } - } + } - for (i = 0; i < ne; i++) - if (mesh -> GetBdrAttribute(i) == mesh -> bdr_attributes[d]) + for (i = 0; i < nelem; i++) + { + mesh->GetBdrElementVertices(elem[i], vertices); + if (cplane == 2) { - mesh->GetBdrElementVertices (i, vertices); - if (cplane == 2) - { - int n = 0; - for (j = 0; j < vertices.Size(); j ++) - if (node_pos[vertices[j]] >= 0.0) - n++; - if (n < vertices.Size()) - continue; // with the next boundary element - } - switch (mesh->GetBdrElementType(i)) - { - case Element::TRIANGLE: - glBegin (GL_TRIANGLES); - break; + int n = 0; + for (j = 0; j < vertices.Size(); j ++) + if (node_pos[vertices[j]] >= 0.0) + n++; + if (n < vertices.Size()) + continue; // with the next boundary element + } + switch (mesh->GetBdrElementType(elem[i])) + { + case Element::TRIANGLE: + glBegin (GL_TRIANGLES); + break; - case Element::QUADRILATERAL: - glBegin (GL_QUADS); - break; - } - mesh->GetBdrPointMatrix (i, pointmat); + case Element::QUADRILATERAL: + glBegin (GL_QUADS); + break; + } + mesh->GetBdrPointMatrix(elem[i], pointmat); - for (j = 0; j < pointmat.Size(); j++) - { - MySetColor ( (*sol)(vertices[j]) , minv , maxv); - glNormal3d (nx(vertices[j]),ny(vertices[j]),nz(vertices[j])); - glVertex3d (pointmat.Elem(0, j), - pointmat.Elem(1, j), - pointmat.Elem(2, j)); - } - glEnd (); + for (j = 0; j < pointmat.Size(); j++) + { + MySetColor((*sol)(vertices[j]), minv ,maxv); + glNormal3d(nx(vertices[j]), ny(vertices[j]), nz(vertices[j])); + glVertex3dv(&pointmat(0, j)); } + glEnd(); + } + } - glEndList (); + glEndList(); } void VisualizationSceneSolution3d::PrepareLines() @@ -1249,7 +1333,8 @@ void VisualizationSceneSolution3d::PrepareLines() Array vertices; - for (i = 0; i < ne; i++){ + for (i = 0; i < ne; i++) + { if (!bdr_attr_to_show[mesh->GetBdrAttribute(i)-1]) continue; mesh -> GetBdrElementVertices (i, vertices); @@ -1282,11 +1367,12 @@ void VisualizationSceneSolution3d::PrepareLines() for (j = 0; j < pointmat.Size(); j++) glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j)); - glEnd (); + glEnd(); break; case 2: - for (j = 0; j < pointmat.Size(); j++) { + for (j = 0; j < pointmat.Size(); j++) + { for (k = 0; k < 3; k++) point[j][k] = pointmat(k,j); point[j][3] = (*sol)(vertices[j]); @@ -1296,7 +1382,7 @@ void VisualizationSceneSolution3d::PrepareLines() } } - glEndList (); + glEndList(); } void VisualizationSceneSolution3d::PrepareLines2() @@ -1308,7 +1394,7 @@ void VisualizationSceneSolution3d::PrepareLines2() glPolygonMode (GL_FRONT_AND_BACK, GL_LINE); int nbe = mesh -> GetNBE(); - DenseMatrix pointmat; + DenseMatrix pointmat, normals; Vector values; RefinedGeometry * RefG; Array vertices; @@ -1340,48 +1426,38 @@ void VisualizationSceneSolution3d::PrepareLines2() // di = GridF -> GetFaceValues (fn, 2, RefG->RefPts, values, pointmat); di = fo % 2; GridF -> GetFaceValues (fn, di, RefG->RefPts, values, pointmat); + ShrinkPoints3D(pointmat, i, fn, fo); - int *RG = &(RefG->RefGeoms[0]); - double pts[][3] = { { pointmat(0,RG[0]), pointmat(1,RG[0]), - pointmat(2,RG[0]) }, - { pointmat(0,RG[1]), pointmat(1,RG[1]), - pointmat(2,RG[1]) }, - { pointmat(0,RG[2]), pointmat(1,RG[2]), - pointmat(2,RG[2]) } }; - double norm[3]; - if (Compute3DUnitNormal (pts[0], pts[1], pts[2], norm)) - cerr << "WARNING: VisualizationSceneSolution3d::PrepareLines2()" - << endl; - if (di != 0 && sc != 0.0) + if (sc != 0.0) { - norm[0] = -norm[0]; - norm[1] = -norm[1]; - norm[2] = -norm[2]; + GetFaceNormals(fn, di, RefG->RefPts, normals); + double norm[3]; + for (int i = 0; i < 3; i++) + norm[i] = 0.0; + for (k = 0; k < normals.Width(); k++) + for (int j = 0; j < 3; j++) + norm[j] += normals(j, k); + double len = sqrt(InnerProd(norm, norm)); + if (len > 0.0) + len = 1.0 / len; + for (int i = 0; i < 3; i++) + norm[i] *= len; + for (int i = 0; i < pointmat.Width(); i++) + { + double val = sc * (values(i) - minv) / (maxv - minv); + for (int j = 0; j < 3; j++) + pointmat(j, i) += val * norm[j]; + } } if (drawmesh == 1) { Array &REdges = RefG->RefEdges; - glBegin (GL_LINES); + glBegin(GL_LINES); for (k = 0; k < REdges.Size(); k++) - { - int *RE = &(REdges[k]); - - if (sc == 0.0) - { - glVertex3d (pointmat(0, RE[0]), pointmat(1, RE[0]), - pointmat(2, RE[0])); - } - else - { - double val = sc * (values(RE[0]) - minv) / (maxv - minv); - glVertex3d (pointmat(0, RE[0])+val*norm[0], - pointmat(1, RE[0])+val*norm[1], - pointmat(2, RE[0])+val*norm[2]); - } - } - glEnd (); + glVertex3dv(&pointmat(0, REdges[k])); + glEnd(); } else if (drawmesh == 2) { @@ -1399,49 +1475,35 @@ void VisualizationSceneSolution3d::PrepareLines2() } for (k = 0; k < RefG->RefGeoms.Size()/sides; k++) { - RG = &(RefG->RefGeoms[k*sides]); + int *RG = &(RefG->RefGeoms[k*sides]); - if (sc == 0.0) - { - for (j = 0; j < sides; j++) - { - for (int ii = 0; ii < 3; ii++) - point[j][ii] = pointmat(ii, RG[j]); - point[j][3] = values(RG[j]); - } - } - else + for (j = 0; j < sides; j++) { - for (j = 0; j < sides; j++) - { - double val = (values(RG[j]) - minv) / (maxv - minv); - val *= sc; - for (int ii = 0; ii < 3; ii++) - point[j][ii] = pointmat(ii, RG[j])+val*norm[ii]; - point[j][3] = values(RG[j]); - } + for (int ii = 0; ii < 3; ii++) + point[j][ii] = pointmat(ii, RG[j]); + point[j][3] = values(RG[j]); } DrawPolygonLevelLines (point[0], sides, level); } } } - glEndList (); + glEndList(); } -void VisualizationSceneSolution3d::CuttingPlaneFunc (int func) +void VisualizationSceneSolution3d::CuttingPlaneFunc(int func) { - int i, j, k, l, m, n = 0; + int i, j, k, m, n, n2; int flag[8], oedges[6]; static const int tet_edges[12]={0,3, 0,2, 0,1, 1,2, 1,3, 2,3}; static const int hex_edges[24] = { 0,1, 1,2, 3,2, 0,3, 4,5, 5,6, 7,6, 4,7, 0,4, 1,5, 2,6, 3,7 }; static const int hex_cutting[24][3] = - { { 3, 5, 7}, { 8, 16, 19}, { 1, 5, 7}, {10, 18, 21}, - { 1, 3, 7}, {12, 20, 23}, { 1, 3, 5}, {14, 17, 22}, - {11, 13, 15}, { 0, 16, 19}, { 9, 13, 15}, { 2, 18, 21}, - { 9, 11, 15}, { 4, 20, 23}, { 9, 11, 13}, { 6, 17, 22}, - { 6, 14, 22}, { 0, 8, 19}, { 0, 8, 16}, { 2, 10, 21}, - { 2, 10, 18}, { 4, 12, 23}, { 4, 12, 20}, { 6, 14, 17} }; + {{ 3, 4, 6}, {17, 9, 18}, { 4, 6, 1}, {19, 11, 20}, + {21, 12, 22}, { 6, 1, 3}, {23, 14, 16}, { 1, 3, 4}, + {18, 0, 17}, {15, 13, 10}, {20, 2, 19}, { 8, 15, 13}, + {10, 8, 15}, {22, 5, 21}, {13, 10, 8}, {16, 7, 23}, + { 9, 18, 0}, { 7, 23, 14}, {11, 20, 2}, { 0, 17, 9}, + {12, 22, 5}, { 2, 19, 11}, {14, 16, 7}, { 5, 21, 12}}; const int *ev; double t, point[6][4], norm[3]; @@ -1450,8 +1512,8 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc (int func) Array nodes; for (i = 0; i < mesh -> GetNE(); i++) { - n = 0; // n will be the number of intersection points - mesh -> GetElementVertices(i,nodes); + n = n2 = 0; // n will be the number of intersection points + mesh -> GetElementVertices(i, nodes); for (j = 0; j < nodes.Size(); j++) if (node_pos[nodes[j]] >= 0.0) flag[j] = 1; @@ -1461,55 +1523,59 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc (int func) { case Element::TETRAHEDRON: ev = tet_edges; - for(j=0; j<6; j++, ev += 2) + for (j = 0; j < 6; j++, ev += 2) if (flag[ev[0]] != flag[ev[1]]) oedges[n++] = j; ev = tet_edges; break; case Element::HEXAHEDRON: + int emark[12]; ev = hex_edges; for (j = 0; j < 12; j++, ev += 2) - if (flag[ev[0]] != flag[ev[1]]) - break; - if (j < 12) + emark[j] = flag[ev[1]] - flag[ev[0]]; + do { - k = 2*j; + for (j = 0; j < 12; j++) + if (emark[j]) + break; + if (j == 12) + break; + k = 2 * j; + if (emark[j] > 0) + k++; do { - for (l = j = 0; j < 3; j++) + for (j = 0; j < 3; j++) { - ev = hex_edges + 2 * (hex_cutting[k][j] / 2); - if (flag[ev[0]] != flag[ev[1]]) - { - m = hex_cutting[k][j]; - l++; + m = hex_cutting[k][j]; + ev = hex_edges + 2 * (m / 2); + if ((m % 2 == 0 && flag[ev[0]] > flag[ev[1]]) || + (m % 2 == 1 && flag[ev[1]] > flag[ev[0]])) break; - } - } - for (j++; j < 3; j++) - { - ev = hex_edges + 2 * (hex_cutting[k][j] / 2); - if (flag[ev[0]] != flag[ev[1]]) - l++; } - if (l != 1 || n == 6) - { - cerr << "VisualizationSceneSolution3d::" - << "CuttingPlaneFunc ("<< func << ")" << endl; - *((int *)NULL) = 0; // force segmentation fault - } - oedges[n++] = k/2; + oedges[n2++] = k/2; + emark[k/2] = 0; k = m; } - while (k/2 != oedges[0]); + while (k/2 != oedges[n]); + if (n == 0) + { + n = n2; + } + else + { + break; + } } + while (1); + n2 -= n; ev = hex_edges; break; default: break; } - if (n > 2) + while (n > 2) { if (shading != 2) mesh -> GetPointMatrix (i, pointmat); @@ -1533,7 +1599,7 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc (int func) node_pos[ nodes[en[1]] ] / ( node_pos[ nodes[en[1]] ] - node_pos[ nodes[en[0]] ] ); - for(k = 0; k < 3; k++) + for (k = 0; k < 3; k++) point[j][k] = t*(pointmat(k,en[0])-pointmat(k,en[1])) + pointmat(k,en[1]); @@ -1545,7 +1611,7 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc (int func) switch (func) { - case 1: // PrepareCuttingPlane() + case 1: // PrepareCuttingPlane() { if (shading == 2) { @@ -1554,20 +1620,35 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc (int func) } else { - if (n > 3) - j = Compute3DUnitNormal (point[0], point[1], point[2], - point[3], norm); - else - j = Compute3DUnitNormal (point[0], point[1], point[2], - norm); + while (1) + { + if (n > 3) + { + j = Compute3DUnitNormal(point[0], point[1], point[2], + point[3], norm); + if (j && n > 4) + { + for (int j = 3; j < n; j++) + for (int i = 0; i < 4; i++) + point[j-2][i] = point[j][i]; + n -= 2; + continue; + } + } + else + j = Compute3DUnitNormal(point[0], point[1], point[2], + norm); + break; + } + if (!j) { - glBegin (GL_POLYGON); + glBegin(GL_POLYGON); glNormal3dv (norm); for(j = 0; j < n; j++) { - MySetColor ( point[j][3] , minv , maxv); - glVertex3dv (point[j]); + MySetColor(point[j][3], minv, maxv); + glVertex3dv(point[j]); } glEnd(); } @@ -1575,7 +1656,7 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc (int func) } break; - case 2: // PrepareCuttingPlaneLines() with mesh + case 2: // PrepareCuttingPlaneLines() with mesh { if (shading == 2) { @@ -1592,7 +1673,7 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc (int func) } break; - case 3: // PrepareCuttingPlaneLines() with level lines + case 3: // PrepareCuttingPlaneLines() with level lines { if (shading == 2) { @@ -1604,12 +1685,17 @@ void VisualizationSceneSolution3d::CuttingPlaneFunc (int func) } break; } + + for (j = 0; j < n2; j++) + oedges[j] = oedges[j+n]; + n = n2; + n2 = 0; } } } -void VisualizationSceneSolution3d::PrepareCuttingPlane(){ - +void VisualizationSceneSolution3d::PrepareCuttingPlane() +{ if (cp_drawelems == 0) return; if (cplane == 0) return; if (cplane == 2) @@ -1621,9 +1707,9 @@ void VisualizationSceneSolution3d::PrepareCuttingPlane(){ glNewList(cplanelist, GL_COMPILE); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); - Set_Material (); + Set_Material(); - CuttingPlaneFunc (1); + CuttingPlaneFunc(1); glEndList(); } @@ -1631,15 +1717,15 @@ void VisualizationSceneSolution3d::PrepareCuttingPlane(){ void VisualizationSceneSolution3d::PrepareCuttingPlane2() { int i, j, n = 0; - double point[4][4], norm[3]; - DenseMatrix pointmat; + double p[4][3], c[4]; + DenseMatrix pointmat, normals; Vector values; RefinedGeometry * RefG; glNewList(cplanelist, GL_COMPILE); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); - Set_Material (); + Set_Material(); double * coord; @@ -1647,9 +1733,9 @@ void VisualizationSceneSolution3d::PrepareCuttingPlane2() Array partition (mesh -> GetNE()); for (i = 0; i < mesh -> GetNE(); i++) { - n = 0; // n will be the number of nodes behind the cutting plane - mesh -> GetElementVertices(i,nodes); - for(j=0; j GetElementVertices(i, nodes); + for (j = 0; j < nodes.Size(); j++) if (node_pos[nodes[j]] >= 0.0) n++; @@ -1667,27 +1753,19 @@ void VisualizationSceneSolution3d::PrepareCuttingPlane2() if (shading != 2) { mesh -> GetFaceVertices (i, nodes); - for(j = 0; j < nodes.Size(); j++) + 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]); + p[j][0] = coord[0]; + p[j][1] = coord[1]; + p[j][2] = coord[2]; + c[j] = (*sol)(nodes[j]); } - if (Compute3DUnitNormal (point[0], point[1], point[2], norm)) - cerr << "WARNING: " - << "VisualizationSceneSolution3d::PrepareCuttingPlane2()" - << endl; - glBegin (GL_POLYGON); - for(j = 0; j < nodes.Size(); j++) - { - MySetColor (point[j][3], minv, maxv); - glNormal3dv ( norm ); - glVertex3dv (point[j]); - } - glEnd(); + if (nodes.Size() == 3) + DrawTriangle(p, c, minv, maxv); + else + DrawQuad(p, c, minv, maxv); } else // shading == 2 { @@ -1697,19 +1775,22 @@ void VisualizationSceneSolution3d::PrepareCuttingPlane2() // and 1 otherwise int di = partition[e1]; GridF -> GetFaceValues (i, di, RefG->RefPts, values, pointmat); + GetFaceNormals(i, di, RefG->RefPts, normals); switch (mesh -> GetFaceBaseGeometry (i)) { case Geometry::TRIANGLE: n = 3; break; case Geometry::SQUARE: n = 4; break; } - DrawRefinedSurf (n, pointmat, values, RefG->RefGeoms); + // DrawRefinedSurf (n, pointmat, values, RefG->RefGeoms); + DrawPatch(pointmat, values, normals, n, RefG->RefGeoms, + minv, maxv, di ? -3 : 2); } // end shading == 2 } glEndList(); } -void VisualizationSceneSolution3d::PrepareCuttingPlaneLines(){ - +void VisualizationSceneSolution3d::PrepareCuttingPlaneLines() +{ if (cp_drawmesh == 0) return; if (cplane == 0) return; if (cplane == 2 && cp_drawmesh != 3) @@ -1730,7 +1811,7 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines(){ case 3: func = 3; break; } - CuttingPlaneFunc (func); + CuttingPlaneFunc(func); glEndList(); } @@ -1752,7 +1833,7 @@ void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2() Array partition (mesh -> GetNE()); for (i = 0; i < mesh -> GetNE(); i++) { - n = 0; // n will be the number of nodes behind the cutting plane + n = 0; // n will be the number of nodes behind the cutting plane mesh -> GetElementVertices(i,nodes); for(j=0; j= 0.0) @@ -1838,12 +1919,12 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() return; } - ne = mesh -> GetNE (); + ne = mesh -> GetNE(); glNewList(lsurflist, GL_COMPILE); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); - Set_Material (); + Set_Material(); levels.SetSize(nlevels); for (l = 0; l < nlevels; l++) @@ -1917,7 +1998,7 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() glVertex3dv (vert[0]); glVertex3dv (vert[1]); glVertex3dv (vert[2]); - glEnd (); + glEnd(); } } else if (j == 2) @@ -1954,7 +2035,7 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() glVertex3dv (vert[1]); glVertex3dv (vert[2]); glVertex3dv (vert[3]); - glEnd (); + glEnd(); } } } @@ -1963,10 +2044,84 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() glEndList(); } -void VisualizationSceneSolution3d :: Draw (){ +void VisualizationSceneSolution3d::ShrinkPoints3D(DenseMatrix &pointmat, + int i, int fn, int fo) +{ + if (shrink != 1.0) + { + int attr = mesh->GetBdrAttribute(i); + for (int k = 0; k < pointmat.Width(); k++) + for (int d = 0; d < 3; d++) + pointmat(d,k) = shrink*pointmat(d,k) + (1-shrink)*bdrc(d,attr-1); + } + + if (shrinkmat != 1.0) + { + int attr, elem1, elem2; + mesh->GetFaceElements(fn, &elem1, &elem2); + if (fo % 2 == 0) + attr = mesh->GetAttribute(elem1); + else + attr = mesh->GetAttribute(elem2); + for (int k = 0; k < pointmat.Width(); k++) + for (int d = 0; d < 3; d++) + pointmat(d,k) = shrinkmat*pointmat(d,k) + (1-shrinkmat)*matc(d,attr-1); + } +} + +void VisualizationSceneSolution3d::ComputeBdrAttrCenter() +{ + DenseMatrix pointmat; + Vector nbdrc(mesh->bdr_attributes.Max()); + + bdrc.SetSize(3,mesh->bdr_attributes.Max()); + bdrc = 0.0; + nbdrc = 0.0; + + for (int i = 0; i < mesh -> GetNBE(); i++) + { + mesh->GetBdrPointMatrix(i, pointmat); + nbdrc(mesh->GetBdrAttribute(i)-1) += pointmat.Width(); + for (int k = 0; k < pointmat.Width(); k++) + for (int d = 0; d < 3; d++) + bdrc(d,mesh->GetBdrAttribute(i)-1) += pointmat(d,k); + } + + for (int i = 0; i < mesh->bdr_attributes.Max(); i++) + if (nbdrc(i) != 0) + for (int d = 0; d < 3; d++) + bdrc(d,i) /= nbdrc(i); +} + +void VisualizationSceneSolution3d::ComputeElemAttrCenter() +{ + DenseMatrix pointmat; + Vector nmatc(mesh->attributes.Max()); + + matc.SetSize(3,mesh->attributes.Max()); + matc = 0.0; + nmatc = 0.0; + + for (int i = 0; i < mesh -> GetNE(); i++) + { + mesh->GetPointMatrix(i, pointmat); + nmatc(mesh->GetAttribute(i)-1) += pointmat.Width(); + for (int k = 0; k < pointmat.Width(); k++) + for (int d = 0; d < 3; d++) + matc(d,mesh->GetAttribute(i)-1) += pointmat(d,k); + } + + for (int i = 0; i < mesh->attributes.Max(); i++) + if (nmatc(i) != 0) + for (int d = 0; d < 3; d++) + matc(d,i) /= nmatc(i); +} + +void VisualizationSceneSolution3d::Draw() +{ glEnable(GL_DEPTH_TEST); - Set_Background (); + Set_Background(); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); // draw colored faces @@ -1998,7 +2153,7 @@ void VisualizationSceneSolution3d :: Draw (){ else DrawColorBar(minv,maxv); - Set_Black_Material (); + Set_Black_Material(); // draw axes if (drawaxes){ glCallList(axeslist); @@ -2021,14 +2176,14 @@ void VisualizationSceneSolution3d :: Draw (){ // tr_eqn[3] = eqn[3] + 1e-2; // glClipPlane(GL_CLIP_PLANE0,tr_eqn); glClipPlane (GL_CLIP_PLANE0, CuttingPlane->Equation()); - Set_Black_Material (); + Set_Black_Material(); glDisable(GL_CLIP_PLANE0); if ( cp_drawmesh ) glCallList(cplanelineslist); glEnable(GL_CLIP_PLANE0); } - Set_Black_Material (); + Set_Black_Material(); glDisable(GL_LIGHTING); // draw lines @@ -2038,7 +2193,7 @@ void VisualizationSceneSolution3d :: Draw (){ glEnable(GL_LIGHTING); if (MatAlpha < 1.0) - Set_Transparency (); + Set_Transparency(); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); if (GetUseTexture()) @@ -2051,7 +2206,7 @@ void VisualizationSceneSolution3d :: Draw (){ { // glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); glCallList (lsurflist); - // Set_Black_Material (); + // Set_Black_Material(); // glPolygonMode (GL_FRONT_AND_BACK, GL_LINE); // glCallList (lsurflist); } @@ -2071,7 +2226,7 @@ void VisualizationSceneSolution3d :: Draw (){ glDisable (GL_TEXTURE_1D); if (MatAlpha < 1.0) - Remove_Transparency (); + Remove_Transparency(); glFlush(); glXSwapBuffers (auxXDisplay(), auxXWindow()); diff --git a/lib/vssolution3d.hpp b/lib/vssolution3d.hpp index e342cee..d4cad93 100644 --- a/lib/vssolution3d.hpp +++ b/lib/vssolution3d.hpp @@ -12,8 +12,6 @@ #ifndef GLVIS_VSSOLUTION_3D #define GLVIS_VSSOLUTION_3D -int Compute3DUnitNormal (double p1[], double p2[], double p3[], double nor[]); - class VisualizationSceneSolution3d : public VisualizationSceneScalarData { protected: @@ -32,6 +30,9 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData void Init(); + void GetFaceNormals(const int FaceNo, const int side, + const IntegrationRule &ir, DenseMatrix &normals); + void DrawRefinedSurf (int n, double *points, int elem, int func, int part = -1); void DrawRefinedSurf (int n, DenseMatrix &pointmat, @@ -48,14 +49,20 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData int TimesToRefine; double FaceShiftScale; - int attr_to_show; Array bdr_attr_to_show; + double shrinkmat; + // Centers of gravity based on the bounday/element attributes + DenseMatrix bdrc, matc; + VisualizationSceneSolution3d(); VisualizationSceneSolution3d(Mesh & m, Vector & s); void SetGridFunction (GridFunction *gf) { GridF = gf; }; + void NewMeshAndSolution(Mesh *new_m, Vector *new_sol, + GridFunction *new_u = NULL, int rescale = 0); + virtual ~VisualizationSceneSolution3d(); virtual void FindNewBox(); @@ -65,13 +72,13 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData virtual void Prepare(); virtual void Draw(); - void ToggleDrawElems () + void ToggleDrawElems() { drawelems = !drawelems; } - void ToggleDrawMesh (); + void ToggleDrawMesh(); - void ToggleShading (); - int GetShading () { return shading; }; + void ToggleShading(); + int GetShading() { return shading; }; virtual void SetShading(int); virtual void SetRefineFactors(int, int); @@ -79,6 +86,7 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData void CuttingPlaneFunc (int type); void CPPrepare(); + void CPMoved(); void PrepareFlat2(); void PrepareLines2(); virtual void PrepareCuttingPlane(); @@ -87,11 +95,11 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData void PrepareCuttingPlaneLines2(); void PrepareLevelSurf(); void ToggleCuttingPlane(); - void ToggleCPDrawElems (); - void ToggleCPDrawMesh (); - void MoveLevelSurf (int); - void NumberOfLevelSurf (int); - virtual void EventUpdateColors () + void ToggleCPDrawElems(); + void ToggleCPDrawMesh(); + void MoveLevelSurf(int); + void NumberOfLevelSurf(int); + virtual void EventUpdateColors() { Prepare(); PrepareCuttingPlane(); PrepareLevelSurf(); if (shading == 2 && drawmesh != 0 && FaceShiftScale != 0.0) @@ -100,6 +108,43 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData virtual void UpdateLevelLines() { PrepareLines(); PrepareCuttingPlaneLines(); } virtual void UpdateValueRange() { } + + /// Shrink the set of points towards attributes center of gravity + void ShrinkPoints3D(DenseMatrix &pointmat, int i, int fn, int fo); + void ComputeBdrAttrCenter(); + void ComputeElemAttrCenter(); }; +inline double InnerProd(double a[], double b[]) +{ + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +inline void CrossProd(double a[], double b[], double cp[]) +{ + cp[0] = a[1] * b[2] - a[2] * b[1]; + cp[1] = a[2] * b[0] - a[0] * b[2]; + cp[2] = a[0] * b[1] - a[1] * b[0]; +} + +inline int Normalize(double v[]) +{ + double len = sqrt(InnerProd(v, v)); + if (len > 0.0) + len = 1.0 / len; + else + return 1; + for (int i = 0; i < 3; i++) + v[i] *= len; + return 0; +} + +int Normalize(DenseMatrix &normals); + +int Compute3DUnitNormal(const double p1[], const double p2[], + const double p3[], double nor[]); + +int Compute3DUnitNormal (const double p1[], const double p2[], + const double p3[], const double p4[], double nor[]); + #endif diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp index 5823183..0b1b5a2 100644 --- a/lib/vsvector.cpp +++ b/lib/vsvector.cpp @@ -11,13 +11,15 @@ #include #include +#include #include #include "mfem.hpp" #include "visual.hpp" -static void VectorKeyHPressed (){ +static void VectorKeyHPressed() +{ cout << endl << "+------------------------------------+" << endl << "| Keys |" << endl @@ -204,24 +206,8 @@ void VisualizationSceneVector::ToggleVectorField() PrepareVectorField(); } -// VisualizationSceneVector::VisualizationSceneVector(istream & in) -// { -// mesh = new DMesh<2>(); -// solx = new Vector(); -// soly = new Vector(); -// -// mesh -> Load (in); -// solx -> Load (in, mesh -> GetNV()); -// soly -> Load (in, mesh -> GetNV()); -// -// sol = new Vector(mesh -> GetNV()); -// -// Init(); -// } - - -VisualizationSceneVector -::VisualizationSceneVector(Mesh & m, Vector & sx, Vector & sy) +VisualizationSceneVector::VisualizationSceneVector(Mesh & m, + Vector & sx, Vector & sy) { mesh = &m; solx = &sx; @@ -234,7 +220,7 @@ ::VisualizationSceneVector(Mesh & m, Vector & sx, Vector & sy) Init(); } -VisualizationSceneVector::VisualizationSceneVector (GridFunction &vgf) +VisualizationSceneVector::VisualizationSceneVector(GridFunction &vgf) { FiniteElementSpace *fes = vgf.FESpace(); if (fes == NULL || fes->GetVDim() != 2) @@ -322,13 +308,13 @@ void VisualizationSceneVector::CycleVec2Scalar() (*sol)(i) = Vec2Scalar((*solx)(i), (*soly)(i)); // update scalar stuff - FindNewBox (); - PrepareAxes (); - PrepareLines (); + FindNewBox(); + PrepareAxes(); + PrepareLines(); DefaultLevelLines(); - PrepareLevelCurves (); + PrepareLevelCurves(); PrepareCP(); - Prepare (); + Prepare(); if (i == 0) maxlen = maxv; @@ -399,7 +385,7 @@ void VisualizationSceneVector::Init() vectorlist = glGenLists(1); displinelist = glGenLists(1); - VisualizationSceneSolution :: Init(); + VisualizationSceneSolution::Init(); PrepareVectorField(); // PrepareDisplacedMesh(); // called by PrepareLines() @@ -430,7 +416,7 @@ void VisualizationSceneVector::Init() maxlen = maxv; } -VisualizationSceneVector::~VisualizationSceneVector () +VisualizationSceneVector::~VisualizationSceneVector() { glDeleteLists (displinelist, 1); glDeleteLists (vectorlist, 1); @@ -528,6 +514,20 @@ void VisualizationSceneVector::GetRefinedValues( } } } + + if (shrink != 1.0) + ShrinkPoints(tr); +} + +int VisualizationSceneVector::GetRefinedValuesAndNormals( + int i, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr, + DenseMatrix &normals) +{ + int have_normals = 0; + + GetRefinedValues(i, ir, vals, tr); + + return have_normals; } void VisualizationSceneVector::PrepareDisplacedMesh() @@ -540,7 +540,7 @@ void VisualizationSceneVector::PrepareDisplacedMesh() // prepare the displaced mesh glNewList(displinelist, GL_COMPILE); - // Set_Material (); + // Set_Material(); // glColor3f (1, 0, 0); glPolygonMode (GL_FRONT_AND_BACK, GL_LINE); @@ -567,7 +567,7 @@ void VisualizationSceneVector::PrepareDisplacedMesh() pointmat.Elem(1, j)+ (*soly)(vertices[j])*(ianim)/ianimmax, zc); - glEnd (); + glEnd(); } } else if (drawdisp < 2) @@ -593,7 +593,7 @@ void VisualizationSceneVector::PrepareDisplacedMesh() glVertex3d (pm(0, RE[k]) + sc * vvals(0, RE[k]), pm(1, RE[k]) + sc * vvals(1, RE[k]), zc); } - glEnd (); + glEnd(); } } else @@ -722,13 +722,13 @@ void VisualizationSceneVector::PrepareDisplacedMesh() } } - glEndList (); + glEndList(); } double new_maxlen; -void VisualizationSceneVector::DrawVector( - double px, double py, double vx, double vy, double cval) +void VisualizationSceneVector::DrawVector(double px, double py, double vx, + double vy, double cval) { double zc = 0.5*(z[0]+z[1]); @@ -828,10 +828,11 @@ void VisualizationSceneVector::PrepareVectorField() while(rerun); } -void VisualizationSceneVector :: Draw (){ +void VisualizationSceneVector::Draw() +{ glEnable(GL_DEPTH_TEST); - Set_Background (); + Set_Background(); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); // model transformation @@ -857,7 +858,7 @@ void VisualizationSceneVector :: Draw (){ else DrawColorBar(minv,maxv); - Set_Black_Material (); + Set_Black_Material(); // draw axes if (drawaxes) @@ -906,7 +907,7 @@ void VisualizationSceneVector :: Draw (){ glCallList(vectorlist); if (MatAlpha < 1.0) - Set_Transparency (); + Set_Transparency(); // draw elements glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); @@ -914,7 +915,7 @@ void VisualizationSceneVector :: Draw (){ glCallList(displlist); if (MatAlpha < 1.0) - Remove_Transparency (); + Remove_Transparency(); if (GetUseTexture()) glDisable (GL_TEXTURE_1D); diff --git a/lib/vsvector.hpp b/lib/vsvector.hpp index 12a0bd4..e172388 100644 --- a/lib/vsvector.hpp +++ b/lib/vsvector.hpp @@ -25,6 +25,9 @@ class VisualizationSceneVector : public VisualizationSceneSolution virtual void GetRefinedValues(int i, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr); + virtual int GetRefinedValuesAndNormals(int i, const IntegrationRule &ir, + Vector &vals, DenseMatrix &tr, + DenseMatrix &normals); double (*Vec2Scalar)(double, double); diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index 94d3445..a026638 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -16,7 +16,8 @@ #include "mfem.hpp" #include "visual.hpp" -static void VectorKeyHPressed (){ +static void VectorKeyHPressed() +{ cout << endl << "+------------------------------------+" << endl << "| Keys |" << endl @@ -60,9 +61,10 @@ static void VectorKeyHPressed (){ << "+------------------------------------+" << endl << "| F1 - X window info and keystrokes |" << endl << "| F2 - Update colors, etc. |" << endl - << "| F3/F4 - Shrink/Zoom subdomains |" << endl + << "| F3/F4 - Shrink/Zoom bdr elements |" << endl << "| F6 - Palete options |" << endl << "| F7 - Manually set min/max value |" << endl + << "| F11/F12 - Shrink/Zoom subdomains |" << endl << "+------------------------------------+" << endl << "| Keypad |" << endl << "+------------------------------------+" << endl @@ -89,14 +91,14 @@ static void VectorKeyHPressed (){ VisualizationSceneVector3d * vsvector3d; extern VisualizationScene * locscene; -//extern VisualizationSceneVector3d * vssol3d; -static void KeyDPressed(){ +static void KeyDPressed() +{ vsvector3d -> ToggleDisplacements(); SendExposeEvent(); } -static void KeyNPressed () +static void KeyNPressed() { if (vsvector3d -> drawdisp) vsvector3d -> ianimd =( (vsvector3d ->ianimd + 1) % @@ -107,7 +109,7 @@ static void KeyNPressed () vsvector3d -> NPressed(); } -static void KeyBPressed () +static void KeyBPressed() { if (vsvector3d -> drawdisp) vsvector3d ->ianimd = ((vsvector3d ->ianimd + @@ -120,9 +122,10 @@ static void KeyBPressed () vsvector3d -> NPressed(); } -static void KeyrPressed (){ +static void KeyrPressed() +{ locscene -> spinning = 0; - auxIdleFunc((void (*)())0); + RemoveIdleFunc(MainLoop); vsvector3d -> CenterObject(); locscene -> ViewAngle = 45.0; locscene -> ViewScale = 1.0; @@ -136,14 +139,15 @@ static void KeyrPressed (){ SendExposeEvent(); } -static void KeyRPressed (){ +static void KeyRPressed() +{ vsvector3d -> ianim = vsvector3d -> ianimd = 0; vsvector3d -> Prepare(); vsvector3d -> PrepareLines(); vsvector3d -> PrepareDisplacedMesh(); } -void VisualizationSceneVector3d::NPressed () +void VisualizationSceneVector3d::NPressed() { if (drawdisp) PrepareDisplacedMesh(); @@ -156,17 +160,19 @@ void VisualizationSceneVector3d::NPressed () SendExposeEvent(); } -static void KeyuPressed(){ +static void KeyuPressed() +{ vsvector3d -> ToggleVectorFieldLevel(+1); SendExposeEvent(); } -static void KeyUPressed(){ +static void KeyUPressed() +{ vsvector3d -> ToggleVectorFieldLevel(-1); SendExposeEvent(); } -void VisualizationSceneVector3d::ToggleVectorFieldLevel (int v) +void VisualizationSceneVector3d::ToggleVectorFieldLevel(int v) { int i; for (i = 0; i < vflevel.Size(); i++) @@ -179,12 +185,14 @@ void VisualizationSceneVector3d::ToggleVectorFieldLevel (int v) vsvector3d -> PrepareVectorField(); } -static void KeywPressed(){ +static void KeywPressed() +{ vsvector3d -> AddVectorFieldLevel(); SendExposeEvent(); } -static void KeyWPressed(){ +static void KeyWPressed() +{ vsvector3d -> RemoveVectorFieldLevel(); SendExposeEvent(); } @@ -205,17 +213,20 @@ void VisualizationSceneVector3d::RemoveVectorFieldLevel() vsvector3d -> PrepareVectorField(); } -static void KeyvPressed(){ +static void KeyvPressed() +{ vsvector3d -> ToggleVectorField(1); SendExposeEvent(); } -static void KeyVPressed(){ +static void KeyVPressed() +{ vsvector3d -> ToggleVectorField(-1); SendExposeEvent(); } -void VisualizationSceneVector3d::ToggleVectorField(int i){ +void VisualizationSceneVector3d::ToggleVectorField(int i) +{ if (drawvector == 0 && i == -1) drawvector = 5; else @@ -223,25 +234,8 @@ void VisualizationSceneVector3d::ToggleVectorField(int i){ PrepareVectorField(); } -// VisualizationSceneVector3d::VisualizationSceneVector3d(istream & in) -// { -// mesh = new DMesh<2>(); -// solx = new Vector(); -// soly = new Vector(); -// solz = new Vector(); -// -// mesh -> Load (in); -// solx -> Load (in, mesh -> GetNV()); -// soly -> Load (in, mesh -> GetNV()); -// solz -> Load (in, mesh -> GetNV()); -// -// sol = new Vector(mesh -> GetNV()); -// -// Init(); -// } - -VisualizationSceneVector3d -::VisualizationSceneVector3d(Mesh & m, Vector & sx, Vector & sy, Vector & sz) +VisualizationSceneVector3d::VisualizationSceneVector3d(Mesh & m, Vector & sx, + Vector & sy, Vector & sz) { mesh = &m; solx = &sx; @@ -256,7 +250,7 @@ ::VisualizationSceneVector3d(Mesh & m, Vector & sx, Vector & sy, Vector & sz) Init(); } -VisualizationSceneVector3d::VisualizationSceneVector3d (GridFunction &vgf) +VisualizationSceneVector3d::VisualizationSceneVector3d(GridFunction &vgf) { FiniteElementSpace *fes = vgf.FESpace(); if (fes == NULL || fes->GetVDim() != 3) @@ -320,7 +314,7 @@ void VisualizationSceneVector3d::Init() vectorlist = glGenLists(1); displinelist = glGenLists(1); - VisualizationSceneSolution3d :: Init(); + VisualizationSceneSolution3d::Init(); PrepareVectorField(); PrepareDisplacedMesh(); @@ -361,7 +355,7 @@ void VisualizationSceneVector3d::Init() } } -VisualizationSceneVector3d::~VisualizationSceneVector3d () +VisualizationSceneVector3d::~VisualizationSceneVector3d() { glDeleteLists (vectorlist, 1); glDeleteLists (displinelist, 1); @@ -379,77 +373,61 @@ VisualizationSceneVector3d::~VisualizationSceneVector3d () } -void VisualizationSceneVector3d :: PrepareFlat (){ +void VisualizationSceneVector3d::PrepareFlat() +{ int i, j; glNewList (displlist, GL_COMPILE); - Set_Material (); + Set_Material(); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); int nbe = mesh -> GetNBE(); DenseMatrix pointmat; Array vertices; + double p[4][3], c[4]; for (i = 0; i < nbe; i++) { if (!bdr_attr_to_show[mesh->GetBdrAttribute(i)-1]) continue; - switch (mesh->GetBdrElementType(i)) - { - case Element::TRIANGLE: - glBegin (GL_TRIANGLES); - break; + mesh->GetBdrPointMatrix(i, pointmat); + mesh->GetBdrElementVertices(i, vertices); - case Element::QUADRILATERAL: - glBegin (GL_QUADS); - break; - } - mesh->GetBdrPointMatrix (i, pointmat); - mesh->GetBdrElementVertices (i, vertices); - - for (j = 0; j < pointmat.Size(); j++) { + for (j = 0; j < pointmat.Width(); j++) + { pointmat(0, j) += (*solx)(vertices[j])*(ianim)/ianimmax; pointmat(1, j) += (*soly)(vertices[j])*(ianim)/ianimmax; pointmat(2, j) += (*solz)(vertices[j])*(ianim)/ianimmax; } - double v10[] = { pointmat(0,1)-pointmat(0,0), - pointmat(1,1)-pointmat(1,0), - pointmat(2,1)-pointmat(2,0) }; - double v21[] = { pointmat(0,2)-pointmat(0,1), - pointmat(1,2)-pointmat(1,1), - pointmat(2,2)-pointmat(2,1) }; - - double norm[] = { v10[1]*v21[2]-v10[2]*v21[1], - v10[2]*v21[0]-v10[0]*v21[2], - v10[0]*v21[1]-v10[1]*v21[0] }; - double rlen = 1.0/sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]); - - glNormal3d (norm[0]*rlen, norm[1]*rlen, norm[2]*rlen); - - for (j = 0; j < pointmat.Size(); j++) { - MySetColor ( (*sol)(vertices[j]) , minv, maxv); - glVertex3d (pointmat(0, j), - pointmat(1, j), - pointmat(2, j)); + for (j = 0; j < pointmat.Width(); j++) + { + p[j][0] = pointmat(0, j); + p[j][1] = pointmat(1, j); + p[j][2] = pointmat(2, j); + c[j] = (*sol)(vertices[j]); } - glEnd (); + if (j == 3) + DrawTriangle(p, c, minv, maxv); + else + DrawQuad(p, c, minv, maxv); } - glEndList (); + glEndList(); } void VisualizationSceneVector3d::PrepareFlat2() { - int i, j, k, fn, fo, ft, di; + int i, j, k, fn, fo, ft, di, have_normals; double bbox_diam, vmin, vmax, mm; int nbe = mesh -> GetNBE(); - DenseMatrix pointmat; + DenseMatrix pointmat, normals; DenseMatrix vec_vals; Vector values; RefinedGeometry * RefG; Array vertices; + double norm[3]; bbox_diam = sqrt ( (x[1]-x[0])*(x[1]-x[0]) + (y[1]-y[0])*(y[1]-y[0]) + @@ -459,7 +437,7 @@ void VisualizationSceneVector3d::PrepareFlat2() glNewList (displlist, GL_COMPILE); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); - Set_Material (); + Set_Material(); ft = 1; for (i = 0; i < nbe; i++) @@ -483,9 +461,19 @@ void VisualizationSceneVector3d::PrepareFlat2() TimesToRefine); // di = GridF -> GetFaceValues (fn, 2, RefG->RefPts, values, pointmat); di = fo % 2; - GridF -> GetFaceValues (fn, di, RefG->RefPts, values, pointmat); - VecGridF -> GetFaceVectorValues (fn, di, RefG->RefPts, - vec_vals, pointmat); + GridF->GetFaceValues(fn, di, RefG->RefPts, values, pointmat); + if (ianim > 0) + { + VecGridF->GetFaceVectorValues(fn, di, RefG->RefPts, vec_vals, + pointmat); + pointmat.Add(double(ianim)/ianimmax, vec_vals); + have_normals = 0; + } + else + { + GetFaceNormals(fn, di, RefG->RefPts, normals); + have_normals = 1; + } if (ft) { vmin = values.Min(); @@ -500,11 +488,10 @@ void VisualizationSceneVector3d::PrepareFlat2() if (mm > vmax) vmax = mm; } - if (ianim > 0) - pointmat.Add (double(ianim)/ianimmax, vec_vals); + ShrinkPoints3D(pointmat, i, fn, fo); int sides; - switch (mesh -> GetBdrElementType (i)) + switch (mesh->GetBdrElementType(i)) { case Element::TRIANGLE: sides = 3; @@ -514,81 +501,38 @@ void VisualizationSceneVector3d::PrepareFlat2() sides = 4; break; } - int *RG = &(RefG->RefGeoms[0]); - double pts[][3] = { { pointmat(0,RG[0]), pointmat(1,RG[0]), - pointmat(2,RG[0]) }, - { pointmat(0,RG[1]), pointmat(1,RG[1]), - pointmat(2,RG[1]) }, - { pointmat(0,RG[2]), pointmat(1,RG[2]), - pointmat(2,RG[2]) } }; - double norm[3]; - Compute3DUnitNormal (pts[0], pts[1], pts[2], norm); - if (di != 0 && sc != 0.0) - { - norm[0] = -norm[0]; - norm[1] = -norm[1]; - norm[2] = -norm[2]; - } - for (k = 0; k < RefG->RefGeoms.Size()/sides; k++) + if (sc != 0.0 && have_normals) { - RG = &(RefG->RefGeoms[k*sides]); - switch (sides) + for (int i = 0; i < 3; i++) + norm[i] = 0.0; + Normalize(normals); + for (k = 0; k < normals.Width(); k++) + for (int j = 0; j < 3; j++) + norm[j] += normals(j, k); + Normalize(norm); + for (int i = 0; i < pointmat.Width(); i++) { - case 3: - glBegin (GL_TRIANGLES); - break; - - case 4: - glBegin (GL_QUADS); - break; + double val = sc * (values(i) - minv) / (maxv - minv); + for (int j = 0; j < 3; j++) + pointmat(j, i) += val * norm[j]; } - - if (sc == 0.0) - { - glNormal3dv (norm); - for (j = 0; j < sides; j++) - { - MySetColor ( values(RG[j]) , minv , maxv); - glVertex3d (pointmat(0, RG[j]), pointmat(1, RG[j]), - pointmat(2, RG[j])); - } - } - else - { - double nnorm[3]; - for (j = 0; j < 3; j++) - { - double val = sc * (values(RG[j]) - minv) / (maxv - minv); - for (int l = 0; l < 3; l++) - pts[j][l] = pointmat(l, RG[j]) + val * norm[l]; - } - Compute3DUnitNormal (pts[0], pts[1], pts[2], nnorm); - glNormal3dv (nnorm); - for (j = 0; j < 3; j++) - { - MySetColor ( values(RG[j]) , minv , maxv); - glVertex3dv (pts[j]); - } - for ( ; j < sides; j++) - { - double val = (values(RG[j]) - minv) / (maxv - minv); - MySetColor ( values(RG[j]) , minv , maxv ); - val *= sc; - glVertex3d (pointmat(0, RG[j])+val*norm[0], - pointmat(1, RG[j])+val*norm[1], - pointmat(2, RG[j])+val*norm[2]); - } - } - glEnd (); + have_normals = 0; } + + have_normals = have_normals ? 2 : 0; + if (di) + have_normals = -1 - have_normals; + DrawPatch(pointmat, values, normals, sides, RefG->RefGeoms, + minv, maxv, have_normals); } - glEndList (); + glEndList(); cout << "VisualizationSceneVector3d::PrepareFlat2() : [min,max] = [" << vmin << "," << vmax << "]" << endl; } -void VisualizationSceneVector3d :: Prepare (){ +void VisualizationSceneVector3d::Prepare() +{ int i,j; switch (shading){ @@ -603,13 +547,14 @@ void VisualizationSceneVector3d :: Prepare (){ } glNewList (displlist, GL_COMPILE); - Set_Material (); + Set_Material(); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); int ne = mesh -> GetNBE(); int nv = mesh -> GetNV(); DenseMatrix pointmat; Array vertices; + double nor[3]; Vector nx(nv); Vector ny(nv); @@ -629,30 +574,26 @@ void VisualizationSceneVector3d :: Prepare (){ mesh->GetBdrPointMatrix (i, pointmat); mesh->GetBdrElementVertices (i, vertices); - for (j = 0; j < pointmat.Size(); j++) { + for (j = 0; j < pointmat.Size(); j++) + { pointmat(0, j) += (*solx)(vertices[j])*(ianim)/ianimmax; pointmat(1, j) += (*soly)(vertices[j])*(ianim)/ianimmax; pointmat(2, j) += (*solz)(vertices[j])*(ianim)/ianimmax; } - double v10[] = { pointmat(0,1)-pointmat(0,0), - pointmat(1,1)-pointmat(1,0), - pointmat(2,1)-pointmat(2,0)}; - double v21[] = { pointmat(0,2)-pointmat(0,1), - pointmat(1,2)-pointmat(1,1), - pointmat(2,2)-pointmat(2,1) }; - double norm[] = { v10[1]*v21[2]-v10[2]*v21[1], - v10[2]*v21[0]-v10[0]*v21[2], - v10[0]*v21[1]-v10[1]*v21[0] }; - double rlen = 1.0/sqrt(norm[0]*norm[0]+norm[1]*norm[1]+ - norm[2]*norm[2]); - - for (j = 0; j < pointmat.Size(); j++) - { - nx(vertices[j]) += norm[0]*rlen; - ny(vertices[j]) += norm[1]*rlen; - nz(vertices[j]) += norm[2]*rlen; - } + if (pointmat.Width() == 3) + j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1), + &pointmat(0,2), nor); + else + j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1), + &pointmat(0,2), &pointmat(0,3), nor); + if (j == 0) + for (j = 0; j < pointmat.Size(); j++) + { + nx(vertices[j]) += nor[0]; + ny(vertices[j]) += nor[1]; + nz(vertices[j]) += nor[2]; + } } for (i = 0; i < ne; i++) @@ -670,25 +611,22 @@ void VisualizationSceneVector3d :: Prepare (){ } mesh->GetBdrPointMatrix (i, pointmat); mesh->GetBdrElementVertices (i, vertices); - - for (j = 0; j < pointmat.Size(); j++) { + for (j = 0; j < pointmat.Size(); j++) + { pointmat(0, j) += (*solx)(vertices[j])*(ianim)/ianimmax; pointmat(1, j) += (*soly)(vertices[j])*(ianim)/ianimmax; pointmat(2, j) += (*solz)(vertices[j])*(ianim)/ianimmax; } - for (j = 0; j < pointmat.Size(); j++) { - MySetColor ( (*sol)(vertices[j]) , minv , maxv); - glNormal3d (nx(vertices[j]),ny(vertices[j]),nz(vertices[j])); - glVertex3d (pointmat(0, j), - pointmat(1, j), - pointmat(2, j)); + MySetColor((*sol)(vertices[j]), minv, maxv); + glNormal3d(nx(vertices[j]), ny(vertices[j]), nz(vertices[j])); + glVertex3dv(&pointmat(0, j)); } - glEnd (); + glEnd(); } } - glEndList (); + glEndList(); } void VisualizationSceneVector3d::PrepareLines() @@ -734,9 +672,9 @@ void VisualizationSceneVector3d::PrepareLines() for (j = 0; j < pointmat.Size(); j++) glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j) ); - glEnd (); + glEnd(); } - glEndList (); + glEndList(); } void VisualizationSceneVector3d::PrepareLines2() @@ -786,6 +724,8 @@ void VisualizationSceneVector3d::PrepareLines2() if (ianim > 0) pointmat.Add (double(ianim)/ianimmax, vec_vals); + ShrinkPoints3D(pointmat, i, fn, fo); + int *RG = &(RefG->RefGeoms[0]); double pts[][3] = { { pointmat(0,RG[0]), pointmat(1,RG[0]), pointmat(2,RG[0]) }, @@ -827,7 +767,7 @@ void VisualizationSceneVector3d::PrepareLines2() } } } - glEnd (); + glEnd(); } else if (drawmesh == 2) { @@ -871,10 +811,11 @@ void VisualizationSceneVector3d::PrepareLines2() } } } - glEndList (); + glEndList(); } -void VisualizationSceneVector3d::PrepareDisplacedMesh(){ +void VisualizationSceneVector3d::PrepareDisplacedMesh() +{ int i, j, ne = mesh -> GetNBE(); DenseMatrix pointmat; Array vertices; @@ -882,7 +823,7 @@ void VisualizationSceneVector3d::PrepareDisplacedMesh(){ // prepare the displaced mesh glNewList(displinelist, GL_COMPILE); - Set_Material (); + Set_Material(); glColor3f (1, 0, 0); glPolygonMode (GL_FRONT_AND_BACK, GL_LINE); @@ -908,12 +849,13 @@ void VisualizationSceneVector3d::PrepareDisplacedMesh(){ for (j = 0; j < pointmat.Size(); j++) glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j) ); - glEnd (); + glEnd(); } - glEndList (); + glEndList(); } -void ArrowsDrawOrNot (Array l[], int nv, Vector & sol, int nl, Array & level) +void ArrowsDrawOrNot (Array l[], int nv, Vector & sol, + int nl, Array & level) { static int first_time = 1; static int nll = nl; @@ -957,8 +899,9 @@ int ArrowDrawOrNot (double v, int nl, Array & level) return 0; } -void VisualizationSceneVector3d::DrawVector (int type, double v0, double v1, double v2, - double sx, double sy, double sz, double s) +void VisualizationSceneVector3d::DrawVector (int type, double v0, double v1, + double v2, double sx, double sy, + double sz, double s) { static int nv = mesh -> GetNV(); static double volume = (x[1]-x[0])*(y[1]-y[0])*(z[1]-z[0]); @@ -1005,13 +948,14 @@ void VisualizationSceneVector3d::DrawVector (int type, double v0, double v1, dou } } -void VisualizationSceneVector3d::PrepareVectorField(){ +void VisualizationSceneVector3d::PrepareVectorField() +{ int i, nv = mesh -> GetNV(); double *vertex; glNewList(vectorlist, GL_COMPILE); - Set_Material (); + Set_Material(); switch (drawvector) { @@ -1059,11 +1003,6 @@ void VisualizationSceneVector3d::PrepareVectorField(){ ArrowsDrawOrNot (l, nv, *sol, nl, level); int j,k; - /* - double volume = (x[1]-x[0])*(y[1]-y[0])*(z[1]-z[0]) - * xscale * yscale * zscale; - double h = pow(volume, 0.333) / 10; - */ for (k = 0; k < vflevel.Size(); k++) { i = vflevel[k]; @@ -1094,12 +1033,12 @@ void VisualizationSceneVector3d::PrepareVectorField(){ } break; } - glEndList (); + glEndList(); } -void VisualizationSceneVector3d::PrepareCuttingPlane(){ - +void VisualizationSceneVector3d::PrepareCuttingPlane() +{ if (cp_drawelems == 0) return; if (cplane == 0) return; if (cplane == 2) @@ -1108,12 +1047,6 @@ void VisualizationSceneVector3d::PrepareCuttingPlane(){ return; } - /* - double volume = (x[1]-x[0])*(y[1]-y[0])*(z[1]-z[0]) - * xscale * yscale * zscale; - double h = pow(volume/nv, 0.333); - */ - int i, j, k, n = 0; int flag[4], ind[6][2]={{0,3},{0,2},{0,1},{1,2},{1,3},{2,3}}; double t, point[4][4], val[4][3]; @@ -1121,7 +1054,7 @@ void VisualizationSceneVector3d::PrepareCuttingPlane(){ glNewList(cplanelist, GL_COMPILE); glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); - Set_Material (); + Set_Material(); DenseMatrix pointmat(3,4); double * coord; @@ -1129,7 +1062,7 @@ void VisualizationSceneVector3d::PrepareCuttingPlane(){ Array nodes; for(i=0; i < mesh -> GetNE(); i++) { - n = 0; // n will be the number of intersection points + n = 0; // n will be the number of intersection points mesh -> GetElementVertices(i,nodes); mesh -> GetPointMatrix (i,pointmat); for(j=0; j<4; j++) @@ -1177,7 +1110,8 @@ void VisualizationSceneVector3d::PrepareCuttingPlane(){ n++; } - if (n > 2){ + if (n > 2) + { double v10[] = { point[1][0]-point[0][0], point[1][1]-point[0][1], point[1][2]-point[0][2] }; @@ -1230,10 +1164,11 @@ void VisualizationSceneVector3d::PrepareCuttingPlane(){ glEndList(); } -void VisualizationSceneVector3d :: Draw (){ +void VisualizationSceneVector3d::Draw() +{ glEnable(GL_DEPTH_TEST); - Set_Background (); + Set_Background(); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); // draw colored faces @@ -1263,7 +1198,7 @@ void VisualizationSceneVector3d :: Draw (){ if (light) glEnable(GL_LIGHTING); - Set_Black_Material (); + Set_Black_Material(); // draw axes if (drawaxes) { @@ -1277,12 +1212,12 @@ void VisualizationSceneVector3d :: Draw (){ { glClipPlane(GL_CLIP_PLANE0,CuttingPlane->Equation()); glEnable(GL_CLIP_PLANE0); - Set_Black_Material (); + Set_Black_Material(); if ( cp_drawmesh ) glCallList(cplanelineslist); } - Set_Black_Material (); + Set_Black_Material(); // draw lines if (drawmesh) @@ -1291,11 +1226,21 @@ void VisualizationSceneVector3d :: Draw (){ if (drawdisp) glCallList(displinelist); + if (drawvector == 1 || drawvector > 3) + glCallList(vectorlist); + if (MatAlpha < 1.0) - Set_Transparency (); + Set_Transparency(); + + if (GetUseTexture()) + { + glEnable (GL_TEXTURE_1D); + glColor4d(1, 1, 1, 1); + } // draw vector field - glCallList(vectorlist); + if (drawvector == 2 || drawvector == 3) + glCallList(vectorlist); // draw elements if (drawelems) @@ -1304,8 +1249,11 @@ void VisualizationSceneVector3d :: Draw (){ if (cplane && cp_drawelems) glCallList(cplanelist); + if (GetUseTexture()) + glDisable (GL_TEXTURE_1D); + if (MatAlpha < 1.0) - Remove_Transparency (); + Remove_Transparency(); glFlush(); glXSwapBuffers (auxXDisplay(), auxXWindow()); diff --git a/lib/vsvector3d.hpp b/lib/vsvector3d.hpp index 7c02f8d..e6e75ba 100644 --- a/lib/vsvector3d.hpp +++ b/lib/vsvector3d.hpp @@ -54,14 +54,14 @@ class VisualizationSceneVector3d : public VisualizationSceneSolution3d virtual void PrepareCuttingPlane(); - void ToggleDisplacements() {drawdisp = (drawdisp+1)%2;}; + void ToggleDisplacements() {drawdisp = (drawdisp+1)%2;}; virtual void Draw(); - virtual void EventUpdateColors () + virtual void EventUpdateColors() { Prepare(); PrepareVectorField(); PrepareCuttingPlane(); }; - void ToggleVectorFieldLevel (int v); + void ToggleVectorFieldLevel(int v); void AddVectorFieldLevel(); void RemoveVectorFieldLevel(); }; diff --git a/makefile b/makefile index 4ea8678..2edd4c2 100644 --- a/makefile +++ b/makefile @@ -13,10 +13,12 @@ CC = g++ OPTS = -O3 DEBUG_OPTS = -g -DGLVIS_DEBUG +DEFINES = -DGLVIS_MULTISAMPLE=4 + # Take screenshots internally with libtiff or externally with xwd? USE_LIBTIFF = NO # Link with LAPACK? (needed if MFEM was compiled with LAPACK support) -USE_LAPACK = NO +USE_LAPACK = YES # GLVis requires the MFEM library MFEM_DIR = ../mfem @@ -49,21 +51,39 @@ TIFF_LIBS_NO = TIFF_OPTS = $(TIFF_OPTS_$(USE_LIBTIFF)) TIFF_LIBS = $(TIFF_LIBS_$(USE_LIBTIFF)) -# Targets - COPTS = $(TIFF_OPTS) $(GL_OPTS) $(OPTS) LIBS = $(MFEM_LIB) $(LAPACK_LIBS) $(X11_LIB) $(GL_LIBS) $(TIFF_LIBS) +CCC = $(CC) $(COPTS) $(DEFINES) + +# generated with 'echo lib/*.c*' +SOURCE_FILES = lib/aux_gl.cpp lib/aux_vis.cpp lib/gl2ps.c lib/material.cpp\ + lib/openglvis.cpp lib/tk.cpp lib/vsdata.cpp lib/vssolution3d.cpp\ + lib/vssolution.cpp lib/vsvector3d.cpp lib/vsvector.cpp +OBJECT_FILES1 = $(SOURCE_FILES:.cpp=.o) +OBJECT_FILES = $(OBJECT_FILES1:.c=.o) +# generated with 'echo lib/*.h*' +HEADER_FILES = lib/aux_gl.hpp lib/aux_vis.hpp lib/gl2ps.h lib/material.hpp\ + lib/openglvis.hpp lib/palettes.hpp lib/tk.h lib/visual.hpp lib/vsdata.hpp\ + lib/vssolution3d.hpp lib/vssolution.hpp lib/vsvector3d.hpp lib/vsvector.hpp + +# Targets + +.SUFFIXES: .c .cpp .o +.cpp.o: + cd $(