diff --git a/CMakeLists.txt b/CMakeLists.txt index 8c56abd9f8f08f40af55407887f0a79b2c193214..e0058faac7f9acb47e7cf30cfe6d7df187c91bcc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -101,6 +101,7 @@ option(COMPILE_TESTS "Compile some tests" OFF) option(USE_PARTITIONER_PARMETIS "Use ParMetis partitioner" OFF) option(USE_PARTITIONER_ZOLTAN "Use Zoltan partitioner" OFF) option(USE_SOLVER_METIS "Use METIS for matrix reordering" OFF) +option(USE_SOLVER_MONDRIAAN "Use Mondriaan for matrix reordering" OFF) option(USE_SOLVER_PETSC "Use PETSc solvers" OFF) option(USE_SOLVER_TRILINOS "Use Trilinos solvers" OFF) #option(USE_AUTODIFF_OPENCL "Use OpenCL for automatic differentiation (under work)" OFF) @@ -174,6 +175,20 @@ if(USE_SOLVER_METIS) endif() endif() +if(USE_SOLVER_MONDRIAAN) + find_package(MONDRIAAN) + if(NOT MONDRIAAN_FOUND) + set(USE_SOLVER_MONDRIAAN OFF) + message("Mondriaan NOT FOUND") + else() + link_directories(${MONDRIAAN_LIBRARY_DIRS}) + include_directories(${MONDRIAAN_INCLUDE_DIRS}) + set(USE_SOLVER_MONDRIAAN ON) + message("Mondriaan FOUND") + endif() +endif() + + if(USE_PARTITIONER_ZOLTAN) find_package(ZOLTAN) if(NOT ZOLTAN_FOUND) diff --git a/Debuggers/VisualStudio/readme.txt b/Debuggers/VisualStudio/readme.txt index cb57ebaffc549aa38a4a075e4dd4da21b0b7a594..50f4a3772f5cbb102bc9e2fa4d7784f591195461 100644 --- a/Debuggers/VisualStudio/readme.txt +++ b/Debuggers/VisualStudio/readme.txt @@ -2,6 +2,6 @@ Copy and paste contents of autoexp.dat into C:\Program Files (x86)\Microsoft Visual Studio 10.0\Common7\Packages\Debugger\autoexp.dat or the directry where Visual Studio was installed. -In future capabilities will grove. +In future capabilities will grow. Tested with Visual studio 2010. diff --git a/cmake_find/LICENSE b/cmake_find/LICENSE new file mode 100644 index 0000000000000000000000000000000000000000..b8bb415782104065ca6eed11a4b17ec04553bedc --- /dev/null +++ b/cmake_find/LICENSE @@ -0,0 +1,23 @@ +Copyright (c) 2014-2015 +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/cmake_find/readme.txt b/cmake_find/readme.txt new file mode 100644 index 0000000000000000000000000000000000000000..236b2dd3ce7489eb89886c8aefb8f49d0e4a7baf --- /dev/null +++ b/cmake_find/readme.txt @@ -0,0 +1,8 @@ +This folder contains scripts for CMake that aid in compilation of the INMOST library with third party libraries. + +The original files were obtained through internet from the following projects: +- PETSc +- MSTK +- FindOpenCL + +And may appear in original or modified form. \ No newline at end of file diff --git a/examples/DrawMatrix/CMakeLists.txt b/examples/DrawMatrix/CMakeLists.txt index ae55446df837038e002d8ab32e129bd8e4254cee..c53110192528849e9e0e3f0682eca387b2e1b129 100644 --- a/examples/DrawMatrix/CMakeLists.txt +++ b/examples/DrawMatrix/CMakeLists.txt @@ -34,6 +34,10 @@ if(OPENGL_FOUND) message("linking DrawMatrix with Metis") target_link_libraries(DrawMatrix ${METIS_LIBRARIES}) endif() + if(USE_SOLVER_MONDRIAAN) + message("linking DrawMatrix with Mondriaan") + target_link_libraries(DrawMatrix ${MONDRIAAN_LIBRARIES}) + endif() endif() diff --git a/examples/FVDiscr/CMakeLists.txt b/examples/FVDiscr/CMakeLists.txt index 06810616a298d978a912d856139ce8f8ca264624..373a1701ef737f3fb915672655e647d75ee13b04 100644 --- a/examples/FVDiscr/CMakeLists.txt +++ b/examples/FVDiscr/CMakeLists.txt @@ -25,6 +25,10 @@ if(USE_SOLVER) message("linking FVDiscr with Metis") target_link_libraries(FVDiscr ${METIS_LIBRARIES}) endif() + if(USE_SOLVER_MONDRIAAN) + message("linking FVDiscr with Mondriaan") + target_link_libraries(FVDiscr ${MONDRIAAN_LIBRARIES}) + endif() endif() diff --git a/examples/FVDiscr/main.cpp b/examples/FVDiscr/main.cpp index e5d94fcb2646bd5e9f277f85732aac3c8a88ddbc..a127844afafb0d0c91a29f9e0e97c9dcca16fad6 100644 --- a/examples/FVDiscr/main.cpp +++ b/examples/FVDiscr/main.cpp @@ -128,6 +128,7 @@ int main(int argc,char ** argv) ttt = Timer(); Solver S(Solver::INNER_ILU2); // Specify the linear solver to ASM+ILU2+BiCGStab one + S.SetParameterReal("absolute_tolerance",1e-8); Solver::Matrix A; // Declare the matrix of the linear system to be solved Solver::Vector x,b; // Declare the solution and the right-hand side vectors @@ -242,6 +243,8 @@ int main(int argc,char ** argv) ttt = Timer(); + Tag error = m->CreateTag("error",DATA_REAL,CELL,NONE,1); + Storage::real err_C = 0.0, err_L2 = 0.0; for( Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell ) if( cell->GetStatus() != Element::Ghost ) @@ -251,6 +254,7 @@ int main(int argc,char ** argv) if (err > err_C) err_C = err; err_L2 += err * err * cell->Volume(); + cell->Real(error) = err; // x[cell->Integer(id)] = err; } err_C = m->AggregateMax(err_C); // Compute the maximal C norm for the error diff --git a/examples/MatSolve/CMakeLists.txt b/examples/MatSolve/CMakeLists.txt index a3f67b6dee7eb04866abae9a66cb02476addb1e6..b1333dfe32d996b653007329410febc1918be5dd 100644 --- a/examples/MatSolve/CMakeLists.txt +++ b/examples/MatSolve/CMakeLists.txt @@ -33,6 +33,10 @@ if(USE_SOLVER) message("linking MatSolve with Metis") target_link_libraries(MatSolve ${METIS_LIBRARIES}) endif() + if(USE_SOLVER_MONDRIAAN) + message("linking MatSolve with Mondriaan") + target_link_libraries(MatSolve ${MONDRIAAN_LIBRARIES}) + endif() endif() diff --git a/examples/OctreeCutcell/octgrid.cpp b/examples/OctreeCutcell/octgrid.cpp index 2f574c79334db9b7238e4905a12a3c10b334a59a..01a12e630a7347139a537dfbd0c38e763f56fc79 100644 --- a/examples/OctreeCutcell/octgrid.cpp +++ b/examples/OctreeCutcell/octgrid.cpp @@ -538,6 +538,7 @@ class incident_matrix } return success; } +/* Storage::real compute_measure(dynarray & data) { Storage::real measure = 0; @@ -622,6 +623,7 @@ class incident_matrix } return measure; } +*/ void recursive_find(unsigned node, unsigned length) { if( !min_loop.empty() && length > min_loop.size() ) return; @@ -851,7 +853,7 @@ public: } } while( min_loop.empty() && first != UINT_MAX ); - for(dynarray::iterator it = min_loop.begin(); it != min_loop.end(); ++it) + for(typename dynarray::iterator it = min_loop.begin(); it != min_loop.end(); ++it) ret.push_back(it->self()); //ret.insert(ret.end(),min_loop.begin(),min_loop.end()); min_loop.clear(); diff --git a/examples/OldDrawGrid/main.cpp b/examples/OldDrawGrid/main.cpp index d9318639677aae09508e23dfc15248f592976d5c..9f627508f4e52b40d864edf3dec6674581d3532c 100644 --- a/examples/OldDrawGrid/main.cpp +++ b/examples/OldDrawGrid/main.cpp @@ -11,6 +11,17 @@ #include #include #include "my_glut.h" +#include + +inline static unsigned int flip(const unsigned int * fp) +{ + unsigned int mask = -((int)(*fp >> 31)) | 0x80000000; + return *fp ^ mask; +} +#define _0(x) (x & 0x7FF) +#define _1(x) (x >> 11 & 0x7FF) +#define _2(x) (x >> 22 ) + using namespace INMOST; @@ -21,7 +32,7 @@ int width = 800, height = 800; double sleft = 1e20, sright = -1e20, sbottom = 1e20, stop = -1e20, sfar = -1e20, snear = 1e20; double shift[3] = {0,0,0}; bool perspective = false; -bool boundary = true, planecontrol = false, clipupdate = false, bndupdate = true, clipboxupdate = false; +bool boundary = true, planecontrol = false, clipupdate = false, bndupdate = true, clipboxupdate = false, draw_volumetric = false; Mesh::GeomParam table; @@ -38,6 +49,13 @@ Mesh::GeomParam table; Storage::real p[3] = {0,0,0}, n[3] = {0,0,1}; ElementArray boundary_faces; +double amplitude = 10; +double radius = 25; +char visualization_prompt[8192]; +bool visualization_prompt_active = false; +Tag visualization_tag; +ElementType visualization_type; + void printtext(const char * fmt, ... ) { @@ -55,18 +73,216 @@ void printtext(const char * fmt, ... ) } } +const int name_width = 32; +const int type_width = 14; +const int elems_width = 10; +const int sparse_width = 10; +const int length_width = 10; + +void PrintTag(Tag t) +{ + std::cout << std::setw(name_width) << t.GetTagName() << std::setw(type_width) << DataTypeName(t.GetDataType()); + int num = 0; + char elems[7] = "NEFCSM"; + std::string print = ""; + for(ElementType etype = NODE; etype <= MESH; etype = etype << 1) + { + if( t.isDefined(etype) ) + { + print = print + elems[ElementNum(etype)]; + num++; + } + } + std::cout << std::setw(elems_width) << print; + print = ""; + for(ElementType etype = NODE; etype <= MESH; etype = etype << 1) + { + if( t.isSparse(etype) ) + { + print = print + elems[ElementNum(etype)]; + num++; + } + } + std::cout << std::setw(sparse_width) << print; + if( t.GetSize() == ENUMUNDEF ) + std::cout << std::setw(length_width) << "VAR" << std::endl; + else + std::cout << std::setw(length_width) << t.GetSize() << std::endl; +} + +void PrintTags(Mesh * m, ElementType etypes) +{ + std::cout << std::setw(name_width) << "Name" << std::setw(type_width) << "Type" << std::setw(elems_width) << "Element" << std::setw(sparse_width) << "Sparse" << std::setw(length_width) << "Length" << std::endl; + for(Mesh::iteratorTag t = m->BeginTag(); t != m->EndTag(); ++t ) + { + bool print = false; + for(ElementType etype = NODE; etype <= MESH; etype = etype << 1) if( (etype&etypes) && t->isDefined(etype) ) {print = true; break;} + if( print ) PrintTag(*t); + } +} + +struct color_t +{ + float c[4]; + color_t() {memset(c,0,sizeof(float)*4);} + color_t(float r, float g, float b) + { + c[0] = r; + c[1] = g; + c[2] = b; + c[3] = 1.0; + } + color_t(float r, float g, float b, float a) + { + c[0] = r; + c[1] = g; + c[2] = b; + c[3] = a; + } + color_t(const color_t & other) {memcpy(c,other.c,sizeof(float)*4);} + color_t & operator =(color_t const & other) + { + memmove(c,other.c,sizeof(float)*4); + return *this; + } + void set_color() const {glColor4fv(c);} + float & r() {return c[0];} + float & g() {return c[1];} + float & b() {return c[2];} + float & a() {return c[3];} + float r() const {return c[0];} + float g() const {return c[1];} + float b() const {return c[2];} + float a() const {return c[3];} + color_t operator *(float mult){return color_t(c[0]*mult,c[1]*mult,c[2]*mult,c[3]*mult);} + color_t operator +(color_t other){return color_t(c[0]+other.c[0],c[1]+other.c[1],c[2]+other.c[2],c[3]+other.c[3]);} +}; + +class color_bar +{ + float min, max; + std::vector ticks; //ticks from 0 to 1 for each color + std::vector colors; //4 floats for each tick + std::string comment; +public: + color_bar() + { + min = 0; + max = 1; + comment = ""; + ticks.push_back(0.f); + ticks.push_back(0.2f); + ticks.push_back(0.4f); + ticks.push_back(0.6f); + ticks.push_back(0.8f); + ticks.push_back(1.f); + colors.push_back(color_t(1,0,0)); + colors.push_back(color_t(1,1,0)); + colors.push_back(color_t(0,1,0)); + colors.push_back(color_t(0,1,1)); + colors.push_back(color_t(0,0,1)); + colors.push_back(color_t(1,0,1)); + + + + + //colors.push_back(color_t(1,0,0)); + } + void set_comment(std::string text) {comment = text;} + void set_min(float newmin) { min = newmin;} + void set_max(float newmax) { max = newmax;} + color_t pick_color(float value) + { + float t = (value-min)/(max-min); + std::vector::iterator it = std::lower_bound(ticks.begin(),ticks.end(),t); + size_t pos = it-ticks.begin(); + if( it == ticks.end() || pos >= ticks.size() ) + { + return colors.back(); + } + if( pos == 0 ) + { + return colors[0]; + } + float interp = (t-ticks[pos-1])/(ticks[pos]-ticks[pos-1]); + return (colors[pos]*interp+colors[pos-1]*(1-interp)); + } + void Draw() + { + float text_pos = -0.89; + float left = -0.95; + float right = -0.9; + float bottom = -0.75; + float top = 0.75; + glBegin(GL_QUADS); + for(int i = 0; i < ticks.size()-1; ++i) + { + colors[i].set_color(); + glVertex2f(left,bottom+ticks[i]*(top-bottom)); + glVertex2f(right,bottom+ticks[i]*(top-bottom)); + colors[i+1].set_color(); + glVertex2f(right,bottom+ticks[(i+1)]*(top-bottom)); + glVertex2f(left,bottom+ticks[(i+1)]*(top-bottom)); + } + glEnd(); + + glColor4f(0,0,0,1); + for(int i = 0; i < ticks.size(); ++i) + { + glRasterPos2f(text_pos,bottom+ticks[i]*(top-bottom)); + printtext("%f",min+ticks[i]*(max-min)); + } + if( comment != "") + { + glRasterPos2f(left,bottom-0.04); + printtext("%s",comment.c_str()); + } + } +} CommonColorBar; + class face2gl { - double dist; + float dist; double c[4]; double cnt[3]; bool flag; ElementType etype; Storage::integer id; std::vector verts; + std::vector colors; + color_t cntcolor; public: - face2gl():verts() {etype = NONE; id = 0; dist = 0; flag = false; memset(c,0,sizeof(double)*4);} + static void radix_sort_dist(std::vector & set) + { + static std::vector tmp; + tmp.resize(set.size()); + unsigned int i; + const unsigned int kHist = 2048; + unsigned int b0[kHist * 3]; + unsigned int *b1 = b0 + kHist; + unsigned int *b2 = b1 + kHist; + memset(b0,0,sizeof(unsigned int)*kHist*3); + for (i = 0; i < set.size(); i++) + { + unsigned int fi = flip((unsigned int *)&set[i].dist); + ++b0[_0(fi)]; ++b1[_1(fi)]; ++b2[_2(fi)]; + } + { + unsigned int sum0 = 0, sum1 = 0, sum2 = 0; + for (i = 0; i < kHist; i++) + { + b0[kHist-1] = b0[i] + sum0; b0[i] = sum0 - 1; sum0 = b0[kHist-1]; + b1[kHist-1] = b1[i] + sum1; b1[i] = sum1 - 1; sum1 = b1[kHist-1]; + b2[kHist-1] = b2[i] + sum2; b2[i] = sum2 - 1; sum2 = b2[kHist-1]; + } + } + for (i = 0; i < set.size(); i++) tmp[++b0[_0(flip((unsigned int *)&set[i].dist))]] = set[i]; + for (i = 0; i < set.size(); i++) set[++b1[_1(flip((unsigned int *)&tmp[i].dist))]] = tmp[i]; + for (i = 0; i < set.size(); i++) tmp[++b2[_2(flip((unsigned int *)&set[i].dist))]] = set[i]; + for (i = 0; i < set.size(); i++) set[i] = tmp[set.size()-1-i]; + } + face2gl():verts(),colors() {etype = NONE; id = 0; dist = 0; flag = false; memset(c,0,sizeof(double)*4);} face2gl(const face2gl & other) :verts(other.verts) { etype = other.etype; @@ -80,6 +296,8 @@ public: c[2] = other.c[2]; c[3] = other.c[3]; flag = other.flag; + colors = other.colors; + cntcolor = other.cntcolor; } face2gl & operator =(face2gl const & other) { @@ -95,17 +313,66 @@ public: c[2] = other.c[2]; c[3] = other.c[3]; flag = other.flag; + colors = other.colors; + cntcolor = other.cntcolor; return *this; } ~face2gl() {} void draw_colour() const { - glColor4dv(c); - for(unsigned k = 0; k < verts.size(); k+=3) + if( colors.empty() ) { - glVertex3dv(cnt); - glVertex3dv(&verts[k]); - glVertex3dv(&verts[(k+3)%verts.size()]); + glColor4dv(c); + for(unsigned k = 0; k < verts.size(); k+=3) + { + glVertex3dv(cnt); + glVertex3dv(&verts[k]); + glVertex3dv(&verts[(k+3)%verts.size()]); + } + } + else + { + for(unsigned k = 0; k < verts.size(); k+=3) + { + cntcolor.set_color(); + glVertex3dv(cnt); + colors[k/3].set_color(); + glVertex3dv(&verts[k]); + colors[(k/3+1)%colors.size()].set_color(); + glVertex3dv(&verts[(k+3)%verts.size()]); + } + } + } + void draw_colour_alpha(double alpha) const + { + if( colors.empty() ) + { + double cc[4] = {c[0],c[1],c[2],alpha}; + glColor4dv(c); + for(unsigned k = 0; k < verts.size(); k+=3) + { + glVertex3dv(cnt); + glVertex3dv(&verts[k]); + glVertex3dv(&verts[(k+3)%verts.size()]); + } + } + else + { + for(unsigned k = 0; k < verts.size(); k+=3) + { + color_t t = cntcolor; + t.a() = alpha; + t.set_color(); + glVertex3dv(cnt); + t = colors[k/3]; + t.a() = alpha; + t.set_color(); + glVertex3dv(&verts[k]); + t = colors[(k/3+1)%colors.size()]; + t.a() = alpha; + t.set_color(); + glVertex3dv(&verts[(k+3)%verts.size()]); + } } } void draw() const @@ -129,13 +396,15 @@ public: void set_color(double r, double g, double b, double a) {c[0] = r; c[1] = g; c[2] = b; c[3] = a;} void add_vert(double x, double y, double z) {unsigned s = (unsigned)verts.size(); verts.resize(s+3); verts[s] = x; verts[s+1] = y; verts[s+2] = z;} void add_vert(double v[3]) {verts.insert(verts.end(),v,v+3);} + void add_color(color_t c) {colors.push_back(c);} double * get_vert(int k) {return &verts[k*3];} unsigned size() {return (unsigned)verts.size()/3;} - void set_center(double _cnt[3]) + void set_center(double _cnt[3], color_t c = color_t(0,0,0,0)) { cnt[0] = _cnt[0]; cnt[1] = _cnt[1]; cnt[2] = _cnt[2]; + cntcolor = c; } void get_center(float _cnt[3]) { @@ -162,6 +431,17 @@ public: cnt[0] /= (verts.size()/3)*1.0; cnt[1] /= (verts.size()/3)*1.0; cnt[2] /= (verts.size()/3)*1.0; + if( !colors.empty() ) compute_center_color(); + } + void compute_center_color() + { + cntcolor.r() = 0; + cntcolor.g() = 0; + cntcolor.b() = 0; + cntcolor.a() = 0; + for(INMOST_DATA_ENUM_TYPE k = 0; k < colors.size(); k++) + cntcolor = cntcolor + colors[k]; + cntcolor = cntcolor*(1.0f/static_cast(colors.size())); } void compute_dist(double cam[3]) { @@ -172,10 +452,26 @@ public: void set_elem(ElementType _etype, Storage::integer _id) {etype = _etype; id = _id;} Element get_elem(Mesh * m) {return m->ElementByLocalID(etype,id);} }; +face2gl DrawFace(Element f); std::vector all_boundary; std::vector clip_boundary; +void draw_faces_nc(std::vector & set, int highlight = -1) +{ + glColor4f(0,1,0,0.1); + glBegin(GL_TRIANGLES); + for(INMOST_DATA_ENUM_TYPE q = 0; q < set.size() ; q++) set[q].draw(); + glEnd(); + if( highlight != -1 ) + { + glColor4f(1,0,0,1); + glBegin(GL_TRIANGLES); + set[highlight].draw(); + glEnd(); + } +} + void draw_faces(std::vector & set, int highlight = -1) { glBegin(GL_TRIANGLES); @@ -190,6 +486,13 @@ void draw_faces(std::vector & set, int highlight = -1) } } +void draw_faces_alpha(std::vector & set, double alpha) +{ + glBegin(GL_TRIANGLES); + for(INMOST_DATA_ENUM_TYPE q = 0; q < set.size() ; q++) set[q].draw_colour_alpha(alpha); + glEnd(); +} + void draw_edges(std::vector & set, int highlight = -1) { glBegin(GL_LINES); @@ -205,6 +508,14 @@ void draw_edges(std::vector & set, int highlight = -1) } +void draw_faces_interactive_nc(std::vector & set) +{ + glColor4f(0,1,0,0.1); + glBegin(GL_TRIANGLES); + for(INMOST_DATA_ENUM_TYPE q = 0; q < set.size() ; q++) if( set[q].get_flag() ) set[q].draw(); + glEnd(); +} + void draw_faces_interactive(std::vector & set) { glBegin(GL_TRIANGLES); @@ -212,6 +523,13 @@ void draw_faces_interactive(std::vector & set) glEnd(); } +void draw_faces_interactive_alpha(std::vector & set, double alpha) +{ + glBegin(GL_TRIANGLES); + for(INMOST_DATA_ENUM_TYPE q = 0; q < set.size() ; q++) if( set[q].get_flag() ) set[q].draw_colour_alpha(alpha); + glEnd(); +} + void draw_edges_interactive(std::vector & set) { glBegin(GL_LINES); @@ -245,6 +563,328 @@ Storage::integer clip_plane_edge(double sp0[3], double sp1[3], double p[3], doub } } +typedef struct point +{ + float coords[3]; + float diam; + float dist; + int id; + /* + point & operator = (point const & b) + { + coords[0] = b.coords[0]; + coords[1] = b.coords[1]; + coords[2] = b.coords[2]; + diam = b.diam; + dist = b.dist; + id = b.id; + } + + point(const point & b) + { + coords[0] = b.coords[0]; + coords[1] = b.coords[1]; + coords[2] = b.coords[2]; + diam = b.diam; + dist = b.dist; + id = b.id; + } + */ +} point_t; + +bool operator <(point_t & a, point_t & b) {return a.dist < b.dist;} + + + + +class volumetric +{ + std::vector points; + Mesh * m; + + void radix_sort_dist(std::vector & set) + { + static std::vector tmp; + tmp.resize(set.size()); + unsigned int i; + const unsigned int kHist = 2048; + unsigned int b0[kHist * 3]; + unsigned int *b1 = b0 + kHist; + unsigned int *b2 = b1 + kHist; + memset(b0,0,sizeof(unsigned int)*kHist*3); + for (i = 0; i < set.size(); i++) + { + unsigned int fi = flip((unsigned int *)&set[i].dist); + ++b0[_0(fi)]; ++b1[_1(fi)]; ++b2[_2(fi)]; + } + { + unsigned int sum0 = 0, sum1 = 0, sum2 = 0; + for (i = 0; i < kHist; i++) + { + b0[kHist-1] = b0[i] + sum0; b0[i] = sum0 - 1; sum0 = b0[kHist-1]; + b1[kHist-1] = b1[i] + sum1; b1[i] = sum1 - 1; sum1 = b1[kHist-1]; + b2[kHist-1] = b2[i] + sum2; b2[i] = sum2 - 1; sum2 = b2[kHist-1]; + } + } + for (i = 0; i < set.size(); i++) tmp[++b0[_0(flip((unsigned int *)&set[i].dist))]] = set[i]; + for (i = 0; i < set.size(); i++) set[++b1[_1(flip((unsigned int *)&tmp[i].dist))]] = tmp[i]; + for (i = 0; i < set.size(); i++) tmp[++b2[_2(flip((unsigned int *)&set[i].dist))]] = set[i]; + for (i = 0; i < set.size(); i++) set[i] = tmp[set.size()-1-i]; + } +public: + + volumetric(Mesh * _m) + { + m = _m; + + points.resize(m->NumberOfCells()); + int q = 0; + for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) + { + Storage::real cnt[3], cntf[3]; + it->Centroid(cnt); + points[q].coords[0] = cnt[0]; + points[q].coords[1] = cnt[1]; + points[q].coords[2] = cnt[2]; + points[q].id = it->LocalID(); + points[q].diam = 0.f; + ElementArray faces = it->getFaces(); + for(ElementArray::iterator f = faces.begin(); f != faces.end(); ++f) + { + f->Centroid(cntf); + Storage::real d = sqrt((cnt[0]-cntf[0])*(cnt[0]-cntf[0])+(cnt[1]-cntf[1])*(cnt[1]-cntf[1])+(cnt[2]-cntf[2])*(cnt[2]-cntf[2])); + if( points[q].diam < d ) points[q].diam = d; + } + ++q; + } + /* + points.reserve(m->NumberOfNodes()*20); + int q = 0; + for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) + { + point_t pnt; + Storage::real_array cnt = it->Coords(); + pnt.coords[0] = cnt[0]; + pnt.coords[1] = cnt[1]; + pnt.coords[2] = cnt[2]; + pnt.id = it->LocalID(); + points.push_back(pnt); + + ElementArray adj = it->getAdjElements(CELL); + for(ElementArray::iterator jt = adj.begin(); jt != adj.end(); ++jt) + { + Storage::real cnt2[3]; + jt->Centroid(cnt2); + const int ncoefs = 1; + const Storage::real coefs[ncoefs] = {0.82}; + for(int j = 0; j < ncoefs; ++j) + { + pnt.coords[0] = cnt[0]*(1-coefs[j]) + cnt2[0]*coefs[j]; + pnt.coords[1] = cnt[1]*(1-coefs[j]) + cnt2[1]*coefs[j]; + pnt.coords[2] = cnt[2]*(1-coefs[j]) + cnt2[2]*coefs[j]; + pnt.id = it->LocalID(); + points.push_back(pnt); + } + } + + } + */ + printf("number of points %d\n",points.size()); + } + void camera(double pos[3], int interactive) + { + if( interactive ) return; + float posf[3]; + posf[0] = pos[0]; + posf[1] = pos[1]; + posf[2] = pos[2]; + for(int k = 0; k < points.size(); ++k) + { + points[k].dist = sqrtf( + (posf[0]-points[k].coords[0])*(posf[0]-points[k].coords[0])+ + (posf[1]-points[k].coords[1])*(posf[1]-points[k].coords[1])+ + (posf[2]-points[k].coords[2])*(posf[2]-points[k].coords[2])); + } + double t = Timer(); + radix_sort_dist(points); + //std::sort(points.rbegin(),points.rend()); + printf("Time to sort %lf\n",Timer()-t); + } + void draw(int interactive) + { + double origin[3], right[3], up[3]; + GLdouble modelview[16],projection[16]; + GLint viewport[4] = {0,0,1,1}; + glGetDoublev(GL_MODELVIEW_MATRIX, modelview); + glGetDoublev(GL_PROJECTION_MATRIX, projection); + GLdouble outx, outy, outz; // Var's to save the answer in + gluUnProject(0.5, 0.5, 0., + modelview, projection, viewport, + &outx, &outy, &outz); + origin[0] = outx; + origin[1] = outy; + origin[2] = outz; + gluUnProject(1.0, 0.5, 0., + modelview, projection, viewport, + &outx, &outy, &outz); + right[0] = outx; + right[1] = outy; + right[2] = outz; + gluUnProject(0.5, 1.0, 0., + modelview, projection, viewport, + &outx, &outy, &outz); + up[0] = outx; + up[1] = outy; + up[2] = outz; + right[0] -= origin[0]; + right[1] -= origin[1]; + right[2] -= origin[2]; + up[0] -= origin[0]; + up[1] -= origin[1]; + up[2] -= origin[2]; + + double l = sqrt(right[0]*right[0]+right[1]*right[1]+right[2]*right[2]); + if( l ) + { + right[0] /= l; + right[1] /= l; + right[2] /= l; + } + l = sqrt(up[0]*up[0]+up[1]*up[1]+up[2]*up[2]); + if( l ) + { + up[0] /= l; + up[1] /= l; + up[2] /= l; + } + + const float alpha = 0.0075f; + const float mult = 1.0f; + const float rmult = 0.7f; + //glPointSize(5.0); + glColor4f(0.5f,0.5f,0.5f,alpha); + glEnable(GL_BLEND); + //glBegin(GL_TRIANGLES); + //glBegin(GL_QUADS); + if( interactive ) + { + for(int k = 0; k < points.size(); ++k) if( k % 100 == 0 ) + { + if( visualization_tag.isValid() ) + { + color_t c = CommonColorBar.pick_color(m->CellByLocalID(points[k].id)->RealDF(visualization_tag)); + c.a() = alpha; + c.set_color(); + } + + glBegin(GL_TRIANGLE_FAN); + glVertex3f(points[k].coords[0],points[k].coords[1],points[k].coords[2]); + glVertex3f(points[k].coords[0]+(right[0]*rmult+up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]*rmult+up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]*rmult+up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(right[0]*rmult*2)*points[k].diam*mult,points[k].coords[1]+(right[1]*rmult*2)*points[k].diam*mult,points[k].coords[2]+(right[2]*rmult*2)*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(right[0]*rmult-up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]*rmult-up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]*rmult-up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]*rmult-up[0])*points[k].diam*mult,points[k].coords[1]+(-right[1]*rmult-up[1])*points[k].diam*mult,points[k].coords[2]+(-right[2]*rmult-up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]*rmult*2)*points[k].diam*mult,points[k].coords[1]+(-right[1]*rmult*2)*points[k].diam*mult,points[k].coords[2]+(-right[2]*rmult*2)*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]*rmult+up[0])*points[k].diam*mult,points[k].coords[1]+(-right[1]*rmult+up[1])*points[k].diam*mult,points[k].coords[2]+(-right[2]*rmult+up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(right[0]*rmult+up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]*rmult+up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]*rmult+up[2])*points[k].diam*mult); + glEnd(); + + /* + glVertex3f(points[k].coords[0]+(right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]+up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-up[0]*2)*points[k].diam*mult,points[k].coords[1]+(-up[1]*2)*points[k].diam*mult,points[k].coords[2]+(-up[2]*2)*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(-right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(-right[2]+up[2])*points[k].diam*mult); + */ + /* + glVertex3f(points[k].coords[0]+(right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]+up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(right[0]-up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]-up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]-up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]-up[0])*points[k].diam*mult,points[k].coords[1]+(-right[1]-up[1])*points[k].diam*mult,points[k].coords[2]+(-right[2]-up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(-right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(-right[2]+up[2])*points[k].diam*mult); + */ + //glVertex3f(points[k].coords[0]+right[0]*points[k].diam*mult,points[k].coords[1]+right[1]*points[k].diam*mult,points[k].coords[2]+right[2]*points[k].diam*mult); + //glVertex3f(points[k].coords[0]+up[0]*points[k].diam*mult,points[k].coords[1]+up[1]*points[k].diam*mult,points[k].coords[2]+up[2]*points[k].diam*mult); + //glVertex3f(points[k].coords[0]-right[0]*points[k].diam*mult,points[k].coords[1]-right[1]*points[k].diam*mult,points[k].coords[2]-right[2]*points[k].diam*mult); + //glVertex3f(points[k].coords[0]-up[0]*points[k].diam*mult,points[k].coords[1]-up[1]*points[k].diam*mult,points[k].coords[2]-up[2]*points[k].diam*mult); + } + } + else + for(int k = 0; k < points.size(); ++k) + { + if( visualization_tag.isValid() ) + { + color_t c = CommonColorBar.pick_color(m->CellByLocalID(points[k].id)->RealDF(visualization_tag)); + c.a() = alpha; + c.set_color(); + } + + glBegin(GL_TRIANGLE_FAN); + glVertex3f(points[k].coords[0],points[k].coords[1],points[k].coords[2]); + glVertex3f(points[k].coords[0]+(right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]+up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(right[0]*2)*points[k].diam*mult,points[k].coords[1]+(right[1]*2)*points[k].diam*mult,points[k].coords[2]+(right[2]*2)*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(right[0]-up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]-up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]-up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]-up[0])*points[k].diam*mult,points[k].coords[1]+(-right[1]-up[1])*points[k].diam*mult,points[k].coords[2]+(-right[2]-up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]*2)*points[k].diam*mult,points[k].coords[1]+(-right[1]*2)*points[k].diam*mult,points[k].coords[2]+(-right[2]*2)*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(-right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(-right[2]+up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]+up[2])*points[k].diam*mult); + glEnd(); + + /* + glVertex3f(points[k].coords[0]+(right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]+up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-up[0]*2)*points[k].diam*mult,points[k].coords[1]+(-up[1]*2)*points[k].diam*mult,points[k].coords[2]+(-up[2]*2)*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(-right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(-right[2]+up[2])*points[k].diam*mult); + */ + /* + glVertex3f(points[k].coords[0]+(right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]+up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(right[0]-up[0])*points[k].diam*mult,points[k].coords[1]+(right[1]-up[1])*points[k].diam*mult,points[k].coords[2]+(right[2]-up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]-up[0])*points[k].diam*mult,points[k].coords[1]+(-right[1]-up[1])*points[k].diam*mult,points[k].coords[2]+(-right[2]-up[2])*points[k].diam*mult); + glVertex3f(points[k].coords[0]+(-right[0]+up[0])*points[k].diam*mult,points[k].coords[1]+(-right[1]+up[1])*points[k].diam*mult,points[k].coords[2]+(-right[2]+up[2])*points[k].diam*mult); + */ + //glVertex3f(points[k].coords[0]+right[0]*points[k].diam*mult,points[k].coords[1]+right[1]*points[k].diam*mult,points[k].coords[2]+right[2]*points[k].diam*mult); + //glVertex3f(points[k].coords[0]+up[0]*points[k].diam*mult,points[k].coords[1]+up[1]*points[k].diam*mult,points[k].coords[2]+up[2]*points[k].diam*mult); + //glVertex3f(points[k].coords[0]-right[0]*points[k].diam*mult,points[k].coords[1]-right[1]*points[k].diam*mult,points[k].coords[2]-right[2]*points[k].diam*mult); + //glVertex3f(points[k].coords[0]-up[0]*points[k].diam*mult,points[k].coords[1]-up[1]*points[k].diam*mult,points[k].coords[2]-up[2]*points[k].diam*mult); + } + //glEnd(); + glDisable(GL_BLEND); + //glPointSize(1.0); + } +} * CommonVolumetricView; + + +/* +class volumetric2 +{ + std::vector faces; + Mesh * m; +public: + volumetric2(Mesh * _m) + { + m = _m; + + faces.reserve(m->NumberOfFaces()); + int q = 0; + INMOST_DATA_ENUM_TYPE pace = std::max(1,std::min(15,(unsigned)m->NumberOfFaces()/100)); + for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) + { + faces.push_back(DrawFace(it->self())); + if( q%pace == 0 ) faces.back().set_flag(true); + q++; + } + printf("number of faces %d\n",faces.size()); + } + void camera(double pos[3], int interactive) + { + if( interactive ) return; + for(int k = 0; k < faces.size(); ++k) + faces[k].compute_dist(pos); + //face2gl::radix_sort_dist(faces); + std::sort(faces.rbegin(),faces.rend()); + } + void draw(int interactive) + { + if( interactive ) draw_faces_interactive_alpha(faces,0.05); + else draw_faces_alpha(faces,0.05); + } +}* CommonVolumetricView; +*/ class kdtree { int marked; @@ -289,14 +929,6 @@ class kdtree float bd = eb->xyz[2]; return (ad > bd) - (ad < bd); } - inline static unsigned int flip(const unsigned int * fp) - { - unsigned int mask = -((int)(*fp >> 31)) | 0x80000000; - return *fp ^ mask; - } -#define _0(x) (x & 0x7FF) -#define _1(x) (x >> 11 & 0x7FF) -#define _2(x) (x >> 22 ) void radix_sort(int dim, struct entry * temp) { unsigned int i; @@ -600,15 +1232,17 @@ class clipper { struct edge_point { + double val; Storage::real xyz[3]; Storage::integer edge; edge_point(){} - edge_point(Storage::real _xyz[3], Storage::integer n) + edge_point(Storage::real _xyz[3], Storage::integer n, float v) { xyz[0] = _xyz[0]; xyz[1] = _xyz[1]; xyz[2] = _xyz[2]; edge = n; + val = v; } bool operator ==(const edge_point& b) const { @@ -621,7 +1255,7 @@ class clipper }; Tag clip_point, clip_state; kdtree * tree; - Tag clips; + Tag clips, clipsv; MarkerType marker; ElementArray cells; Mesh * mm; @@ -632,6 +1266,7 @@ public: for(INMOST_DATA_ENUM_TYPE k = 0; k < cells.size(); k++) cells[k]->RemMarker(marker); mm->ReleaseMarker(marker); mm->DeleteTag(clips); + mm->DeleteTag(clipsv); mm->DeleteTag(clip_point); mm->DeleteTag(clip_state); } @@ -642,9 +1277,24 @@ public: tree = new kdtree(m); marker = m->CreateMarker(); clips = m->CreateTag("CLIPS",DATA_REAL,CELL,CELL); + clipsv = m->CreateTag("CLIPS_VAL",DATA_REAL,CELL,CELL); clip_point = m->CreateTag("CLIP_POINT",DATA_REAL,EDGE,NONE,3); clip_state = m->CreateTag("CLIP_STATE",DATA_INTEGER,EDGE,NONE,1); } + double compute_value(Edge e, Storage::real * pnt) + { + if( visualization_tag.isValid() ) + { + Storage::real_array c1 = e->getBeg()->Coords(); + Storage::real_array c2 = e->getEnd()->Coords(); + Storage::real d1,d2,t; + d1 = sqrt((pnt[0]-c1[0])*(pnt[0]-c1[0])+(pnt[1]-c1[1])*(pnt[1]-c1[1])+(pnt[2]-c1[2])*(pnt[2]-c1[2])); + d2 = sqrt((c2[0]-c1[0])*(c2[0]-c1[0])+(c2[1]-c1[1])*(c2[1]-c1[1])+(c2[2]-c1[2])*(c2[2]-c1[2])); + t = d1/d2; //(pnt == c2, t = 1 : pnt == c1, t = 0) + return e->getBeg()->RealDF(visualization_tag)*(1-t)+e->getEnd()->RealDF(visualization_tag)*t; + } + else return 0.f; + } void clip_plane(Storage::real p[3], Storage::real n[3]) { const bool print = false; @@ -653,6 +1303,7 @@ public: if( cells[k]->GetMarker(marker) ) { cells[k]->RealArray(clips).clear(); + cells[k]->RealArray(clipsv).clear(); cells[k]->RemMarker(marker); } tree->intersect_plane_edge(clip_point,clip_state,cells,marker,p,n); @@ -679,8 +1330,8 @@ public: if( state == CLIP_FULL ) { nfulledges++; - edge_point n1 = edge_point(&edges[r].getBeg()->Coords()[0],ntotedges); - edge_point n2 = edge_point(&edges[r].getEnd()->Coords()[0],ntotedges); + edge_point n1 = edge_point(&edges[r].getBeg()->Coords()[0],ntotedges,visualization_tag.isValid() ? edges[r].getBeg()->RealDF(visualization_tag) : 0.f); + edge_point n2 = edge_point(&edges[r].getEnd()->Coords()[0],ntotedges,visualization_tag.isValid() ? edges[r].getEnd()->RealDF(visualization_tag) : 0.f); if( npoints % 2 == 0 ) //all privious edges are closed, just add this one { clipcoords.push_back(n1); @@ -707,7 +1358,7 @@ public: } else if( state == CLIP_ENDP ) { - edge_point n = edge_point(&edges[r].RealArrayDF(clip_point)[0],ntotedges); + edge_point n = edge_point(&edges[r].RealArrayDF(clip_point)[0],ntotedges,compute_value(edges[r],&edges[r].RealArrayDF(clip_point)[0])); bool add = true; if( last_edge_type == CLIP_ENDP ) { @@ -743,7 +1394,7 @@ public: } else if( state == CLIP_NODE ) { - edge_point n = edge_point(&edges[r].RealArrayDF(clip_point)[0],ntotedges); + edge_point n = edge_point(&edges[r].RealArrayDF(clip_point)[0],ntotedges,compute_value(edges[r],&edges[r].RealArrayDF(clip_point)[0])); if( print ) { printf("added: "); @@ -783,13 +1434,16 @@ public: { ElementArray nodes = full_face->getNodes(); Storage::real_array cl = cells[k]->RealArray(clips); + Storage::real_array clv = cells[k]->RealArray(clipsv); cl.resize(static_cast(3*nodes.size())); + clv.resize(static_cast(nodes.size())); for(INMOST_DATA_ENUM_TYPE r = 0; r < nodes.size(); r++) { Storage::real_array p = nodes[r].Coords(); cl[0+3*r] = p[0]; cl[1+3*r] = p[1]; cl[2+3*r] = p[2]; + clv[r] = visualization_tag.isValid() ? nodes[r].RealDF(visualization_tag) : 0.0; } cells[k]->SetMarker(marker); } @@ -831,12 +1485,15 @@ public: __FILE__,__LINE__,ntotedges,loopcoords.size()); } Storage::real_array cl = cells[k]->RealArray(clips); + Storage::real_array clv = cells[k]->RealArray(clipsv); cl.resize(static_cast(3*loopcoords.size())); + clv.resize(static_cast(loopcoords.size())); for(INMOST_DATA_ENUM_TYPE r = 0; r < loopcoords.size(); ++r) { cl[r*3+0] = loopcoords[r].xyz[0]; cl[r*3+1] = loopcoords[r].xyz[1]; cl[r*3+2] = loopcoords[r].xyz[2]; + clv[r] = loopcoords[r].val; } loopcoords.clear(); @@ -854,7 +1511,12 @@ public: face2gl f; f.set_color(0.6,0.6,0.6,1); Storage::real_array cl = cells[k]->RealArray(clips); + Storage::real_array clv = cells[k]->RealArray(clipsv); for(INMOST_DATA_ENUM_TYPE q = 0; q < cl.size(); q+=3) f.add_vert(&cl[q]); + if( visualization_tag.isValid() ) + { + for(INMOST_DATA_ENUM_TYPE q = 0; q < clv.size(); q++) f.add_color(CommonColorBar.pick_color(clv[q])); + } f.compute_center(); f.set_elem(cells[k]->GetElementType(),cells[k]->LocalID()); if( k%pace == 0 ) f.set_flag(true); @@ -867,20 +1529,27 @@ public: for(INMOST_DATA_ENUM_TYPE k = 0; k < cells.size(); k+=pace) if( cells[k]->GetMarker(marker)) { Storage::real_array cl = cells[k]->RealArray(clips); + Storage::real_array clv = cells[k]->RealArray(clipsv); Storage::real cnt[3] = {0,0,0}; + Storage::real cntv = 0.0; for(INMOST_DATA_ENUM_TYPE q = 0; q < cl.size(); q+=3) { cnt[0] += cl[q+0]; cnt[1] += cl[q+1]; cnt[2] += cl[q+2]; + cntv += clv[q/3]; } - cnt[0] /= (cl.size()/3); - cnt[1] /= (cl.size()/3); - cnt[2] /= (cl.size()/3); + cnt[0] /= static_cast(cl.size()/3); + cnt[1] /= static_cast(cl.size()/3); + cnt[2] /= static_cast(cl.size()/3); + cntv /= static_cast(cl.size()/3); for(INMOST_DATA_ENUM_TYPE q = 0; q < cl.size(); q+=3) { + if( visualization_tag.isValid() ) CommonColorBar.pick_color(cntv).set_color(); glVertex3dv(cnt); + if( visualization_tag.isValid() ) CommonColorBar.pick_color(clv[q/3]).set_color(); glVertex3dv(&cl[q]); + if( visualization_tag.isValid() ) CommonColorBar.pick_color(clv[(q/3+1)%clv.size()]).set_color(); glVertex3dv(&cl[(q+3)%cl.size()]); } } @@ -927,6 +1596,18 @@ public: for(INMOST_DATA_ENUM_TYPE k = 0; k < nfaces; k++) faces[k] = _faces[k]; tree = new kdtree(mm,faces,nfaces); } + double compute_value(Node n1, Node n2, Storage::real * c1, Storage::real * c2, Storage::real * pnt) + { + if( visualization_tag.isValid() ) + { + Storage::real d1,d2,t; + d1 = sqrt((pnt[0]-c1[0])*(pnt[0]-c1[0])+(pnt[1]-c1[1])*(pnt[1]-c1[1])+(pnt[2]-c1[2])*(pnt[2]-c1[2])); + d2 = sqrt((c2[0]-c1[0])*(c2[0]-c1[0])+(c2[1]-c1[1])*(c2[1]-c1[1])+(c2[2]-c1[2])*(c2[2]-c1[2])); + t = d1/d2; //(pnt == c2, t = 1 : pnt == c1, t = 0) + return n1->RealDF(visualization_tag)*(1-t)+n2->RealDF(visualization_tag)*t; + } + else return 0.f; + } void clip_plane(Storage::real p[3], Storage::real n[3]) { tree->intersect_plane_face(clip_state,p,n); @@ -942,7 +1623,11 @@ public: ElementArray nodes = Face(mm,faces[k])->getNodes(); face2gl f; f.set_color(0.6,0.6,0.6,1); - for(INMOST_DATA_ENUM_TYPE q = 0; q < nodes.size(); q++) f.add_vert(&nodes[q].Coords()[0]); + for(INMOST_DATA_ENUM_TYPE q = 0; q < nodes.size(); q++) + { + f.add_vert(&nodes[q].Coords()[0]); + if( visualization_tag.isValid() ) f.add_color(CommonColorBar.pick_color(nodes[q].RealDF(visualization_tag))); + } f.compute_center(); f.set_elem(GetHandleElementType(faces[k]),GetHandleID(faces[k])); if( k%pace == 0 ) f.set_flag(true); @@ -968,13 +1653,18 @@ public: { coords = nodes[q].Coords(); f.add_vert(&coords[0]); + if( visualization_tag.isValid() ) f.add_color(CommonColorBar.pick_color(nodes[q].RealDF(visualization_tag))); } if( nodepos[q] != nodepos[(q+1)%nodes.size()] ) { Storage::real_array sp0 = nodes[q].Coords(); Storage::real_array sp1 = nodes[(q+1)%nodes.size()].Coords(); Storage::real node[3]; - if( clip_plane_edge(&sp0[0],&sp1[0],p,n,node) > CLIP_NONE) f.add_vert(node); + if( clip_plane_edge(&sp0[0],&sp1[0],p,n,node) > CLIP_NONE) + { + f.add_vert(node); + if( visualization_tag.isValid() ) f.add_color(CommonColorBar.pick_color(compute_value(nodes[q],nodes[(q+1)%nodes.size()],&sp0[0],&sp1[0],node))); + } } } f.compute_center(); @@ -1007,7 +1697,12 @@ public: { ElementArray nodes = Face(mm,faces[k])->getNodes(); glBegin(GL_POLYGON); - for(INMOST_DATA_ENUM_TYPE q = 0; q < nodes.size(); q++) glVertex3dv(&nodes[q].Coords()[0]); + for(INMOST_DATA_ENUM_TYPE q = 0; q < nodes.size(); q++) + { + if( visualization_tag.isValid() ) CommonColorBar.pick_color(nodes[q].RealDF(visualization_tag)).set_color(); + glVertex3dv(&nodes[q].Coords()[0]); + + } glEnd(); } mm->IntegerDF(faces[k],clip_state) = state; @@ -1283,8 +1978,8 @@ void set_matrix3d() else { const double pi = 3.1415926535897932384626433832795; - const double znear = 0.01; - const double zfar = 10000.0; + const double znear = 1; + const double zfar = 1000000.0; const double fH = znear * tan( 60.0 / 360.0 * pi); double fW = fH * aspect; glFrustum(-fW,fW,-fH,fH,znear,zfar); @@ -1378,7 +2073,10 @@ void myclickmotion(int nmx, int nmy) // Mouse clickmotion(nmx,nmy); if( planecontrol ) { - rotatevector((double *)n); + //quatget(n); + reverse_rotatevector_from_stack((double *)n); + reverse_rotatevector((double *)n); + rotatevector_from_stack((double *)n); clipupdate = true; quatinit(); } @@ -1441,11 +2139,97 @@ void myclick(int b, int s, int nmx, int nmy) // Mouse } +class Input +{ +public: + enum InputType {Double,Integer,String}; +private: + std::string str; + std::string comment; + void * input_link; + InputType type; + bool done; + bool canceled; + +public: + Input(int * val, std::string comment) : comment(comment) {input_link = val; type = Integer; canceled = false; done = false; str = "";} + Input(double * val, std::string comment) : comment(comment) {input_link = val; type = Double; canceled = false; done = false; str = "";} + Input(char * val, std::string comment) : comment(comment) {input_link = val; type = String; canceled = false; done = false; str = "";} + Input(void * link, InputType type, std::string comment) : comment(comment), input_link(link), type(type) {canceled = false; done = false; str = "";} + Input(const Input & other):comment(comment), input_link(other.input_link), str(other.str), type(other.type), canceled(other.canceled), done(other.done) {} + Input & operator =(Input const & other) {comment = other.comment; input_link = other.input_link; str = other.str; type = other.type; canceled = other.canceled; done = other.done; return *this;} + ~Input() {} + void KeyPress(char c) + { + if( c == 13 ) + { + + done = true; + if( type == Double ) *((double *)input_link) = atof(str.c_str()); + else if( type == Integer ) *((int *)input_link) = atoi(str.c_str()); + else if( type == String ) strcpy((char *)input_link,str.c_str()); + glutPostRedisplay(); + } + else if( c == 8 ) + { + if( !str.empty() ) str.erase(str.size()-1); + glutPostRedisplay(); + } + else if( c == 27 ) + { + canceled = true; + done = true; + glutPostRedisplay(); + } + else if( type == String || ( (c >= '0' && c <= '9') || ((str.empty() || tolower(*str.rbegin()) == 'e') && c=='+' || c=='-') || (type == Double && (c=='.' || c=='e' || c == 'E'))) ) + { + str += c; + glutPostRedisplay(); + } + } + bool Done() {return done;} + bool Canceled() {return canceled;} + void Draw() + { + float h = 24.0f/(float)height; + + glColor3f(1,1,1); + glBegin(GL_QUADS); + glVertex3f(-0.99,-0.99,1); + glVertex3f(-0.99,-0.99+h,1); + glVertex3f(0.99,-0.99+h,1); + glVertex3f(0.99,-0.99,1); + glEnd(); + glColor3f(0,0,0); + glBegin(GL_LINE_LOOP); + glVertex3f(-0.99,-0.99,1); + glVertex3f(-0.99,-0.99+h,1); + glVertex3f(0.99,-0.99+h,1); + glVertex3f(0.99,-0.99,1); + glEnd(); + + glColor4f(0,0,0,1); + glRasterPos2d(-0.985,-0.985); + //printtext(str.c_str()); + char oldval[4096]; + if( type == Double ) sprintf(oldval,"%g",*(double*)input_link); + else if( type == Integer ) sprintf(oldval,"%g",*(int*)input_link); + else if( type == String ) sprintf(oldval,"%s",(char *)input_link); + printtext("input number (%s[%s]:%s): %s",comment.c_str(),oldval,type == Integer ? "integer": (type == Double ? "double" : "string"), str.c_str()); + } +} * CommonInput = NULL; + + void keyboard(unsigned char key, int x, int y) { (void) x; (void) y; + if( CommonInput != NULL ) + { + CommonInput->KeyPress(key); + return; + } if( key == 27 ) { if( oclipper ) delete oclipper; @@ -1574,6 +2358,11 @@ void keyboard(unsigned char key, int x, int y) glutPostRedisplay(); } + else if( key == 'c' ) + { + draw_volumetric = !draw_volumetric; + glutPostRedisplay(); + } else if( key == 'b' ) { boundary = !boundary; @@ -1589,6 +2378,107 @@ void keyboard(unsigned char key, int x, int y) } else quatpop(); } + else if( key == 'n' ) + { + if( CommonInput == NULL ) CommonInput = new Input(&litude, "Amplitude"); + glutPostRedisplay(); + } + else if( key == 'm' ) + { + if( CommonInput == NULL ) CommonInput = new Input(&radius, "Radius"); + glutPostRedisplay(); + } + else if( key == ',' || key == '<' ) + { + std::cout << "positive shift" << std::endl; + Storage::real radius2 = radius*radius; + for(Mesh::iteratorNode it = mesh->BeginNode(); it != mesh->EndNode(); ++it) + { + Storage::real_array coords = it->Coords(); + Storage::real proj[3] = {0,0,0}; + Storage::real d = -(n[0]*(coords[0]-p[0]) + n[1]*(coords[1]-p[1]) + n[2]*(coords[2]-p[2])); + proj[0] = coords[0] + d*n[0]; + proj[1] = coords[1] + d*n[1]; + proj[2] = coords[2] + d*n[2]; + d = (n[0]*(proj[0]-p[0]) + n[1]*(proj[1]-p[1]) + n[2]*(proj[2]-p[2])); + assert(fabs(d) < 1.0e-7); // check that it is indeed a projection + Storage::real dist = sqrt((p[0]-proj[0])*(p[0]-proj[0])+(p[1]-proj[1])*(p[1]-proj[1])+(p[2]-proj[2])*(p[2]-proj[2])); + Storage::real delta = exp(-dist*dist/radius2); + coords[0] += n[0]*delta*amplitude; + coords[1] += n[1]*delta*amplitude; + coords[2] += n[2]*delta*amplitude; + } + all_boundary.clear(); + INMOST_DATA_ENUM_TYPE pace = std::max(1,std::min(15,(unsigned)boundary_faces.size()/100)); + for(INMOST_DATA_ENUM_TYPE k = 0; k < boundary_faces.size(); k++) + { + all_boundary.push_back(DrawFace(boundary_faces[k])); + if( k%pace == 0 ) all_boundary.back().set_flag(true); + } + if( oclipper ) delete oclipper; + if( bclipper ) delete bclipper; + bclipper = new bnd_clipper(mesh,boundary_faces.data(),(int)boundary_faces.size()); + oclipper = new clipper(mesh); + clipupdate = true; + glutPostRedisplay(); + } + else if( key == '.' || key == '>' ) + { + std::cout << "negative shift" << std::endl; + Storage::real radius2 = radius*radius; + for(Mesh::iteratorNode it = mesh->BeginNode(); it != mesh->EndNode(); ++it) + { + Storage::real_array coords = it->Coords(); + Storage::real proj[3] = {0,0,0}; + Storage::real d = -(n[0]*(coords[0]-p[0]) + n[1]*(coords[1]-p[1]) + n[2]*(coords[2]-p[2])); + proj[0] = coords[0] + d*n[0]; + proj[1] = coords[1] + d*n[1]; + proj[2] = coords[2] + d*n[2]; + d = (n[0]*(proj[0]-p[0]) + n[1]*(proj[1]-p[1]) + n[2]*(proj[2]-p[2])); + assert(fabs(d) < 1.0e-7); // check that it is indeed a projection + Storage::real dist = sqrt((p[0]-proj[0])*(p[0]-proj[0])+(p[1]-proj[1])*(p[1]-proj[1])+(p[2]-proj[2])*(p[2]-proj[2])); + Storage::real delta = exp(-dist*dist/radius2); + coords[0] -= n[0]*delta*amplitude; + coords[1] -= n[1]*delta*amplitude; + coords[2] -= n[2]*delta*amplitude; + } + all_boundary.clear(); + INMOST_DATA_ENUM_TYPE pace = std::max(1,std::min(15,(unsigned)boundary_faces.size()/100)); + for(INMOST_DATA_ENUM_TYPE k = 0; k < boundary_faces.size(); k++) + { + all_boundary.push_back(DrawFace(boundary_faces[k])); + if( k%pace == 0 ) all_boundary.back().set_flag(true); + } + if( oclipper ) delete oclipper; + if( bclipper ) delete bclipper; + bclipper = new bnd_clipper(mesh,boundary_faces.data(),(unsigned)boundary_faces.size()); + oclipper = new clipper(mesh); + clipupdate = true; + glutPostRedisplay(); + } + else if( key == 'p' ) + { + perspective = !perspective; + glutPostRedisplay(); + } + else if( key == 'v' ) + { + if( CommonInput == NULL ) + { + PrintTags(mesh,CELL|FACE|EDGE|NODE); + + CommonInput = new Input(visualization_prompt, "Enter data for visualization as Element:Name:Component"); + visualization_prompt_active = true; + clipupdate = true; + if( visualization_tag.isValid() ) visualization_tag = mesh->DeleteTag(visualization_tag); + } + glutPostRedisplay(); + } + else if( key == 'q' ) + { + mesh->Save("mesh.vtk"); + mesh->Save("mesh.pmf"); + } } void keyboard2(unsigned char key, int x, int y) @@ -1597,6 +2487,7 @@ void keyboard2(unsigned char key, int x, int y) (void) y; if( key == '=' || key == '+' || key == '_' || key == '-' || key == 'w' || key == 's' || key == 'a' || key == 'd' || key == 'r' || key == 'p' || key == 'z') { + //printf("depressed\n"); interactive = false; glutPostRedisplay(); } @@ -1734,6 +2625,7 @@ void draw() whereami(campos[0],campos[1],campos[2]); int picked = -1; + //glTranslated((l+r)*0.5,(b+t)*0.5,(near+far)*0.5); { @@ -1758,7 +2650,7 @@ void draw() clipboxupdate = false; if( current_picker != NULL ) {delete current_picker; current_picker = NULL;} - current_picker = new picker(clip_boundary); + if( !clip_boundary.empty() ) current_picker = new picker(clip_boundary); } @@ -1798,6 +2690,15 @@ void draw() } } } + + + if( draw_volumetric ) + { + if( bndupdate ) CommonVolumetricView->camera(campos,interactive); + CommonVolumetricView->draw(interactive); + } + + if( boundary ) { glEnable(GL_BLEND); @@ -1805,19 +2706,21 @@ void draw() { for(INMOST_DATA_ENUM_TYPE q = 0; q < all_boundary.size() ; q++) all_boundary[q].compute_dist(campos); - std::sort(all_boundary.rbegin(),all_boundary.rend()); - bndupdate = false; + //std::sort(all_boundary.rbegin(),all_boundary.rend()); + face2gl::radix_sort_dist(all_boundary); } glColor4f(0,0,0,0.25); if( interactive ) draw_edges_interactive(all_boundary); else draw_edges(all_boundary); - if( interactive ) draw_faces_interactive(all_boundary); - else draw_faces(all_boundary); + if( interactive ) draw_faces_interactive_nc(all_boundary); + else draw_faces_nc(all_boundary); glDisable(GL_BLEND); } + + if( !interactive && bndupdate) bndupdate = false; } if( picked != -1 ) @@ -1827,7 +2730,7 @@ void draw() set_matrix2d(); - double top = 0.96, left = 0.11, interval = 0.04, bottom = top; + double top = 0.96, left = 0.25, interval = 0.04, bottom = top; Element e = clip_boundary[picked].get_elem(mesh); for(Mesh::iteratorTag t = mesh->BeginTag(); t != mesh->EndTag(); ++t) if( t->isDefined(e->GetElementType()) ) @@ -1869,7 +2772,7 @@ void draw() Storage::integer_array arr = e->IntegerArray(*t); for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++) { - sprintf(temp,"%d %s",arr[k],str); + sprintf(temp,"%s %d",str,arr[k]); strcpy(str,temp); } break; @@ -1879,7 +2782,7 @@ void draw() Storage::real_array arr = e->RealArray(*t); for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++) { - sprintf(temp,"%lf %s",arr[k],str); + sprintf(temp,"%s %lf",str,arr[k]); strcpy(str,temp); } break; @@ -1889,7 +2792,7 @@ void draw() Storage::bulk_array arr = e->BulkArray(*t); for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++) { - sprintf(temp,"%d %s",arr[k],str); + sprintf(temp,"%s %d",str,arr[k]); strcpy(str,temp); } break; @@ -1899,8 +2802,8 @@ void draw() Storage::reference_array arr = e->ReferenceArray(*t); for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++) { - if(arr.at(k) == InvalidHandle()) sprintf(temp,"NULL %s",str); - else sprintf(temp,"%s:%d %s",ElementTypeName(arr[k]->GetElementType()),arr[k]->LocalID(),str); + if(arr.at(k) == InvalidHandle()) sprintf(temp,"%s NULL",str); + else sprintf(temp,"%s %s:%d",str,ElementTypeName(arr[k]->GetElementType()),arr[k]->LocalID()); strcpy(str,temp); } break; @@ -1913,12 +2816,167 @@ void draw() top -= interval; } } + glEnable(GL_DEPTH_TEST); + } + + if( CommonInput != NULL ) + { + glDisable(GL_DEPTH_TEST); + glLoadIdentity(); + set_matrix2d(); + CommonInput->Draw(); + if( CommonInput->Done() ) + { + if( ! CommonInput->Canceled() ) + { + if( visualization_prompt_active ) + { + char typen[1024],name[1024]; + unsigned comp; + int k = 0,l, slen = (int)strlen(visualization_prompt); + for(k = 0; k < slen; ++k) + { + if( visualization_prompt[k] == ':' ) + { + visualization_prompt[k] = '\0'; + break; + } + } + for(l = k+1; l < slen; ++l) + { + if( visualization_prompt[l] == ':' ) + { + visualization_prompt[l] = '\0'; + break; + } + } + + if( k < slen && l < slen && l+1 < slen ) + { + + strcpy(typen,visualization_prompt); + strcpy(name,visualization_prompt+k+1); + comp = atoi(visualization_prompt+l+1); + visualization_prompt[k] = ':'; + visualization_prompt[l] = ':'; + printf("type %s name %s comp %d\n",typen,name,comp); + std::string stype(typen), sname(name); + if( mesh->HaveTag(sname) ) + { + Tag source_tag = mesh->GetTag(sname); + if( source_tag.GetDataType() == DATA_REAL || source_tag.GetDataType() == DATA_INTEGER ) + { + if( comp >= 0 && comp < source_tag.GetSize() ) + { + visualization_type = NONE; + for(size_t q = 0; q < stype.size(); ++q) + { + stype[q] = tolower(stype[q]); + typen[q] = tolower(typen[q]); + } + if( stype == "node" ) visualization_type = NODE; + else if ( stype == "edge" ) visualization_type = EDGE; + else if ( stype == "face" ) visualization_type = FACE; + else if ( stype == "cell" ) visualization_type = CELL; + + if( visualization_type != NONE ) + { + if( source_tag.isDefined(visualization_type) ) + { + float min = 1.0e20, max = -1.0e20; + printf("prepearing data for visualization\n"); + visualization_tag = mesh->CreateTag("VISUALIZATION_TAG",DATA_REAL,NODE,NONE,1); + + for(Mesh::iteratorNode it = mesh->BeginNode(); it != mesh->EndNode(); ++it) + { + ElementArray elems = it->getAdjElements(visualization_type); + Storage::real_array coords = it->Coords(); + Storage::real cnt[3], dist, wgt; + Storage::real val = 0.0, vol = 0.0, res; + for(ElementArray::iterator jt = elems.begin(); jt != elems.end(); ++jt) if( jt->HaveData(source_tag) && jt->GetDataSize(source_tag) > comp ) + { + jt->Centroid(cnt); + dist = (cnt[0]-coords[0])*(cnt[0]-coords[0])+(cnt[1]-coords[1])*(cnt[1]-coords[1])+(cnt[2]-coords[2])*(cnt[2]-coords[2]); + wgt = 1.0/(dist+1.0e-8); + val += wgt * (source_tag.GetDataType() == DATA_REAL ? jt->RealArray(source_tag)[comp] : jt->IntegerArray(source_tag)[comp] ); + vol += wgt; + } + res = val/vol; + if( res < min ) min = res; + if( res > max ) max = res; + it->RealDF(visualization_tag) = res; + } + visualization_tag = mesh->CreateTag("VISUALIZATION_TAG",DATA_REAL,CELL,NONE,1); + for(Mesh::iteratorCell it = mesh->BeginCell(); it != mesh->EndCell(); ++it) + { + ElementArray elems = it->getAdjElements(visualization_type); + Storage::real coords[3]; + it->Centroid(coords); + Storage::real cnt[3], dist, wgt; + Storage::real val = 0.0, vol = 0.0, res; + for(ElementArray::iterator jt = elems.begin(); jt != elems.end(); ++jt) if( jt->HaveData(source_tag) && jt->GetDataSize(source_tag) > comp ) + { + jt->Centroid(cnt); + dist = (cnt[0]-coords[0])*(cnt[0]-coords[0])+(cnt[1]-coords[1])*(cnt[1]-coords[1])+(cnt[2]-coords[2])*(cnt[2]-coords[2]); + wgt = 1.0/(dist+1.0e-8); + val += wgt * (source_tag.GetDataType() == DATA_REAL ? jt->RealArray(source_tag)[comp] : jt->IntegerArray(source_tag)[comp] ); + vol += wgt; + } + res = val/vol; + if( res < min ) min = res; + if( res > max ) max = res; + it->RealDF(visualization_tag) = res; + } + + CommonColorBar.set_min(min); + CommonColorBar.set_max(max); + char comment[1024]; + sprintf(comment,"%s[%d] on %s",name,comp,typen); + CommonColorBar.set_comment(comment); + clipupdate = true; + /* + all_boundary.clear(); + INMOST_DATA_ENUM_TYPE pace = std::max(1,std::min(15,(unsigned)boundary_faces.size()/100)); + for(INMOST_DATA_ENUM_TYPE k = 0; k < boundary_faces.size(); k++) + { + all_boundary.push_back(DrawFace(boundary_faces[k])); + if( k%pace == 0 ) all_boundary.back().set_flag(true); + } + */ + } + else printf("tag %s is not defined on element type %s\n",name, typen); + } + else printf("do not understand element type %s, should be: node, edge, face, cell\n",typen); + } + else printf("component is out of range for tag %s of size %u\n",name,source_tag.GetSize()); + } + else printf("tag %s is not real or integer\n",name); + } + else printf("mesh do not have tag with name %s\n",name); + } + else printf("malformed string %s for visualization\n",visualization_prompt); + visualization_prompt_active = false; + //visualization_prompt[0] = '\0'; + + } + glutPostRedisplay(); + } + delete CommonInput; + CommonInput = NULL; + } + glEnable(GL_DEPTH_TEST); + } + if( visualization_tag.isValid() ) + { + glDisable(GL_DEPTH_TEST); + glLoadIdentity(); + set_matrix2d(); + CommonColorBar.Draw(); glEnable(GL_DEPTH_TEST); } - glutSwapBuffers(); } @@ -1934,6 +2992,9 @@ int main(int argc, char ** argv) mesh = new Mesh(); printf("Started loading mesh.\n"); tt = Timer(); + mesh->SetCommunicator(INMOST_MPI_COMM_WORLD); + mesh->SetParallelFileStrategy(0); + mesh->SetParallelStrategy(1); mesh->SetFileOption("VERBOSITY","2"); if( argc < 2 ) { @@ -1959,6 +3020,15 @@ int main(int argc, char ** argv) //printf("Delete %lg\n",Timer()-tt); //return 0; + std::map elems; + + for(Mesh::iteratorElement it = mesh->BeginElement(CELL|FACE|EDGE|NODE); it != mesh->EndElement(); ++it) + elems[it->GetGeometricType()]++; + + for(std::map::iterator it = elems.begin(); it != elems.end(); ++it ) + std::cout << Element::GeometricTypeName(it->first) << ": " << it->second << std::endl; + + for(ElementType mask = NODE; mask <= CELL; mask = mask << 1) std::cout << ElementTypeName(mask) << " " << mesh->NumberOf(mask) << std::endl; @@ -2001,6 +3071,12 @@ int main(int argc, char ** argv) } printf("Done. Time %lg\n",Timer()-tt); + if( boundary_faces.empty() ) + { + printf("Haven't found any boundary elements of the mesh. Nothing to display.\n"); + return -1; + } + printf("Prepearing set of boundary faces for drawing.\n"); tt = Timer(); INMOST_DATA_ENUM_TYPE pace = std::max(1,std::min(15,(unsigned)boundary_faces.size()/100)); @@ -2011,6 +3087,7 @@ int main(int argc, char ** argv) } printf("Done. Time %g\n",Timer() - tt); + printf("Prepearing interactive mesh clipper.\n"); //tt = Timer(); oclipper = new clipper(mesh); @@ -2018,7 +3095,8 @@ int main(int argc, char ** argv) clipupdate = true; //printf("Done. Time %g\n",Timer() - tt); - + CommonVolumetricView = new volumetric(mesh); + //CommonVolumetricView = new volumetric2(mesh); shift[0] = -(sleft+sright)*0.5; shift[1] = -(sbottom+stop)*0.5; @@ -2051,6 +3129,7 @@ int main(int argc, char ** argv) glutReshapeFunc(reshape); glutKeyboardFunc(keyboard); + glutKeyboardUpFunc(keyboard2); glutMouseFunc(myclick); glutMotionFunc(myclickmotion); glutPassiveMotionFunc(mymotion); diff --git a/examples/OldDrawGrid/rotate.cpp b/examples/OldDrawGrid/rotate.cpp index 0ec534a376adf40dc05e5094c4d13f59f0fffbae..bd259ffc2fb9e6123b1a36e9fc88f6a86efd12c5 100644 --- a/examples/OldDrawGrid/rotate.cpp +++ b/examples/OldDrawGrid/rotate.cpp @@ -168,6 +168,46 @@ void rotatevector(double * vec) vec[2] = ret[2]/ret[3]; } +void reverse_rotatevector(double * vec) +{ + int i; + double rot[16]; + double temp[4] = {vec[0],vec[1],vec[2],1.0}; + double ret[4]; + q.x = -q.x; + q.y = -q.y; + q.z = -q.z; + rot[ 0] = (q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z); + rot[ 1] = 2.*(q.x*q.y - q.w*q.z); + rot[ 2] = 2.*(q.x*q.z + q.w*q.y); + rot[ 3] = 0.0; + rot[ 4] = 2.*(q.x*q.y + q.w*q.z); + rot[ 5] = (q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z); + rot[ 6] = 2.*(q.y*q.z - q.w*q.x); + rot[ 7] = 0.0; + rot[ 8] = 2.*(q.x*q.z - q.w*q.y); + rot[ 9] = 2.*(q.y*q.z + q.w*q.x); + rot[10] = (q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z); + rot[11] = 0.0; + rot[12] = 0.0; + rot[13] = 0.0; + rot[14] = 0.0; + rot[15] = (q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z); + for(i=0; i < 4; i++) + { + ret[i] = temp[0] * rot[i*4]; + ret[i] += temp[1] * rot[i*4+1]; + ret[i] += temp[2] * rot[i*4+2]; + ret[i] += temp[3] * rot[i*4+3]; + } + vec[0] = ret[0]/ret[3]; + vec[1] = ret[1]/ret[3]; + vec[2] = ret[2]/ret[3]; + q.x = -q.x; + q.y = -q.y; + q.z = -q.z; +} + void revrotatevector(double * vec) { int i; @@ -240,6 +280,51 @@ void rotatevector_from_stack(double * vec) vec[2] = ret[2]/ret[3]; } +void quatget(double *vec) +{ + vec[0] = q.x / q.w; + vec[1] = q.y / q.w; + vec[2] = q.z / q.w; +} + +void reverse_rotatevector_from_stack(double * vec) +{ + int i; + struct quaternion q = storage.back(); + q.x = -q.x; + q.y = -q.y; + q.z = -q.z; + double rot[16]; + double temp[4] = {vec[0],vec[1],vec[2],1.0}; + double ret[4]; + rot[ 0] = (q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z); + rot[ 1] = 2.*(q.x*q.y - q.w*q.z); + rot[ 2] = 2.*(q.x*q.z + q.w*q.y); + rot[ 3] = 0.0; + rot[ 4] = 2.*(q.x*q.y + q.w*q.z); + rot[ 5] = (q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z); + rot[ 6] = 2.*(q.y*q.z - q.w*q.x); + rot[ 7] = 0.0; + rot[ 8] = 2.*(q.x*q.z - q.w*q.y); + rot[ 9] = 2.*(q.y*q.z + q.w*q.x); + rot[10] = (q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z); + rot[11] = 0.0; + rot[12] = 0.0; + rot[13] = 0.0; + rot[14] = 0.0; + rot[15] = (q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z); + for(i=0; i < 4; i++) + { + ret[i] = temp[0] * rot[i*4]; + ret[i] += temp[1] * rot[i*4+1]; + ret[i] += temp[2] * rot[i*4+2]; + ret[i] += temp[3] * rot[i*4+3]; + } + vec[0] = ret[0]/ret[3]; + vec[1] = ret[1]/ret[3]; + vec[2] = ret[2]/ret[3]; +} + void rotate() { double rot[16]; diff --git a/examples/OldDrawGrid/rotate.h b/examples/OldDrawGrid/rotate.h index 17977ce28cbe511fd7afd529ef685208ad64d53f..4926701e3a14aa2317b6b4f7a636e70d536ff856 100644 --- a/examples/OldDrawGrid/rotate.h +++ b/examples/OldDrawGrid/rotate.h @@ -20,13 +20,16 @@ void clickmotion(int nmx, int nmy); void motion(int nmx, int nmy); void click(int b, int s, int nmx, int nmy); +void quatget(double * vec); void quatinit(); void quatpush(); void quatpop(); void rotate(); void rotate_from_stack(); void rotatevector(double * vec); +void reverse_rotatevector(double * vec); void revrotatevector(double * vec); void rotatevector_from_stack(double * vec); +void reverse_rotatevector_from_stack(double * vec); #endif diff --git a/examples/Solver/CMakeLists.txt b/examples/Solver/CMakeLists.txt index 0012eec0ca7a03152dbe641a982053a2a0f09264..4463e40575effe99c555254f35878b5dbf47084d 100644 --- a/examples/Solver/CMakeLists.txt +++ b/examples/Solver/CMakeLists.txt @@ -25,6 +25,10 @@ if(USE_SOLVER) message("linking Solver with Metis") target_link_libraries(Solver ${METIS_LIBRARIES}) endif() + if(USE_SOLVER_MONDRIAAN) + message("linking Solver with Mondriaan") + target_link_libraries(Solver ${MONDRIAAN_LIBRARIES}) + endif() endif(USE_SOLVER) if(USE_PARTITIONER) diff --git a/geometry.cpp b/geometry.cpp index 8d9df354c4284a06570f427bfc76872277f2ca67..9f70ff8055ecdb0567dedb28c6c50110508cc84d 100644 --- a/geometry.cpp +++ b/geometry.cpp @@ -529,6 +529,9 @@ namespace INMOST void Mesh::GetGeometricData(HandleType e, GeometricData type, Storage::real * ret) { + assert(e != InvalidHandle()); + assert(ret != NULL); + assert(type == MEASURE || type == CENTROID || type == BARYCENTER || type == NORMAL); ElementType etype = GetHandleElementType(e); integer edim = Element::GetGeometricDimension(GetGeometricType(e)); integer mdim = GetDimensions(); diff --git a/inmost-config.cmake.in b/inmost-config.cmake.in index 46c53737d3746d98faad6664b0b4c014c4c10da5..b6527e886e1108b17434fed6ff1a869bfec93d5d 100644 --- a/inmost-config.cmake.in +++ b/inmost-config.cmake.in @@ -27,6 +27,8 @@ set(INMOST_LIBRARIES inmost) set(USE_MPI @USE_MPI@) set(USE_PARTITIONER_ZOLTAN @USE_PARTITIONER_ZOLTAN@) set(USE_PARTITIONER_PARMETIS @USE_PARTITIONER_PARMETIS@) +set(USE_SOLVER_MONDRIAAN @USE_SOLVER_MONDRIAAN@) +set(USE_SOLVER_METIS @USE_SOLVER_METIS@) set(USE_SOLVER_TRILINOS @USE_SOLVER_TRILINOS@) set(USE_SOLVER_PETSC @USE_SOLVER_PETSC@) @@ -49,6 +51,16 @@ if( USE_PARTITIONER_PARMETIS ) list(APPEND INMOST_INCLUDE_DIRS "@METIS_INCLUDE_DIR@") endif( USE_PARTITIONER_PARMETIS ) +if( USE_SOLVER_MONDRIAAN ) + list(APPEND INMOST_LIBRARIES "@MONDRIAAN_LIBRARIES@") + list(APPEND INMOST_INCLUDE_DIRS "@MONDRIAAN_INCLUDE_DIRS@") + list(APPEND INMOST_LIBRARY_DIRS "@MONDRIAAN_LIBRARY_DIRS@") +endif( USE_SOLVER_MONDRIAAN ) + +if( USE_SOLVER_METIS ) + list(APPEND INMOST_LIBRARIES "@METIS_LIBRARIES@") + list(APPEND INMOST_INCLUDE_DIRS "@METIS_INCLUDE_DIR@") +endif( USE_SOLVER_METIS ) if( USE_SOLVER_TRILINOS ) list(APPEND INMOST_LIBRARIES "@Trilinos_LIBRARIES@") diff --git a/inmost_mesh.h b/inmost_mesh.h index d408bfa2a5e5912f1ffada86efd884e9539a73c8..492cdb5352f995dde4056c4e9c8e7454ab19bde3 100644 --- a/inmost_mesh.h +++ b/inmost_mesh.h @@ -137,7 +137,7 @@ namespace INMOST __INLINE INMOST_DATA_INTEGER_TYPE GetHandleElementNum (HandleType h) {return h >> handle_etype_shift;} __INLINE ElementType GetHandleElementType (HandleType h) {return 1 << GetHandleElementNum(h);} __INLINE HandleType ComposeHandle (ElementType etype, INMOST_DATA_INTEGER_TYPE ID) {return ID == -1 ? InvalidHandle() : ((ElementNum(etype) << handle_etype_shift) + (1+ID));} - __INLINE HandleType ComposeHandle (INMOST_DATA_INTEGER_TYPE etypenum, INMOST_DATA_INTEGER_TYPE ID) {return ID == -1 ? InvalidHandle() : ((etypenum << handle_etype_shift) + (1+ID));} + __INLINE HandleType ComposeHandleNum (INMOST_DATA_INTEGER_TYPE etypenum, INMOST_DATA_INTEGER_TYPE ID) {return ID == -1 ? InvalidHandle() : ((etypenum << handle_etype_shift) + (1+ID));} __INLINE bool isValidHandle (HandleType h) {return h != 0;} @@ -2356,7 +2356,7 @@ namespace INMOST void Exit (); int & GetFuncID () {return func_id;} std::fstream & GetStream (); - std::fstream & WriteTab (std::fstream & f); + std::ostream & WriteTab (std::ostream & f); void FinalizeFile (); static void AtExit (void); #endif @@ -2933,7 +2933,7 @@ namespace INMOST void EndSequentialCode (); //iterator.cpp:::::::::::::::::::::::::::::::::::::::::::::::::: public: - Element ElementByLocalID (integer etypenum, integer lid) {assert(etypenum < 5 && (lid >= 0 && lid < static_cast(links[etypenum].size())) || (etypenum == 5 && lid == 0)); return Element(this,ComposeHandle(etypenum,lid));} + Element ElementByLocalID (integer etypenum, integer lid) {assert(etypenum < 5 && (lid >= 0 && lid < static_cast(links[etypenum].size())) || (etypenum == 5 && lid == 0)); return Element(this,ComposeHandleNum(etypenum,lid));} Element ElementByLocalID (ElementType etype, integer lid) {return ElementByLocalID(ElementNum(etype),lid);} Element ElementByHandle (HandleType h) {return Element(this,h);} @@ -2941,16 +2941,16 @@ namespace INMOST HandleType PrevHandle (HandleType h) const; //returns InvalidHandle() when go beyond first element HandleType NextHandle (HandleType h, ElementType mask) const; HandleType PrevHandle (HandleType h, ElementType mask) const; //returns InvalidHandle() when go beyond first element - HandleType FirstHandle () const {return ComposeHandle(ElementNum(NODE),0);} - HandleType LastHandle () const {return ComposeHandle(ElementNum(MESH),1);} - HandleType FirstHandle (ElementType etype) const {return ComposeHandle(ElementNum(etype),0);} - HandleType LastHandle (ElementType etype) const {integer num = ElementNum(etype); return ComposeHandle(num,static_cast(links[num].size()));} + HandleType FirstHandle () const {return ComposeHandleNum(ElementNum(NODE),0);} + HandleType LastHandle () const {return ComposeHandleNum(ElementNum(MESH),1);} + HandleType FirstHandle (ElementType etype) const {return ComposeHandleNum(ElementNum(etype),0);} + HandleType LastHandle (ElementType etype) const {integer num = ElementNum(etype); return ComposeHandleNum(num,static_cast(links[num].size()));} - Node NodeByLocalID (integer lid) { assert(lid >= 0 && lid < static_cast(links[0].size())); return Node(this,ComposeHandle(0,lid)); } - Edge EdgeByLocalID (integer lid) { assert(lid >= 0 && lid < static_cast(links[1].size())); return Edge(this,ComposeHandle(1,lid)); } - Face FaceByLocalID (integer lid) { assert(lid >= 0 && lid < static_cast(links[2].size())); return Face(this,ComposeHandle(2,lid));} - Cell CellByLocalID (integer lid) { assert(lid >= 0 && lid < static_cast(links[3].size())); return Cell(this,ComposeHandle(3,lid)); } - ElementSet EsetByLocalID (integer lid) { assert(lid >= 0 && lid < static_cast(links[4].size())); return ElementSet(this,ComposeHandle(4,lid)); } + Node NodeByLocalID (integer lid) { assert(lid >= 0 && lid < static_cast(links[0].size())); return Node(this,ComposeHandleNum(0,lid)); } + Edge EdgeByLocalID (integer lid) { assert(lid >= 0 && lid < static_cast(links[1].size())); return Edge(this,ComposeHandleNum(1,lid)); } + Face FaceByLocalID (integer lid) { assert(lid >= 0 && lid < static_cast(links[2].size())); return Face(this,ComposeHandleNum(2,lid));} + Cell CellByLocalID (integer lid) { assert(lid >= 0 && lid < static_cast(links[3].size())); return Cell(this,ComposeHandleNum(3,lid)); } + ElementSet EsetByLocalID (integer lid) { assert(lid >= 0 && lid < static_cast(links[4].size())); return ElementSet(this,ComposeHandleNum(4,lid)); } integer NodeNextLocalID (integer lid) const {++lid; while(lid < static_cast(links[0].size()) && links[0][lid] == -1) ++lid; return lid;} integer EdgeNextLocalID (integer lid) const {++lid; while(lid < static_cast(links[1].size()) && links[1][lid] == -1) ++lid; return lid;} diff --git a/inmost_options_cmake.h b/inmost_options_cmake.h index 2c6b35a07a74341e4270322efda3eb531ed1258d..76165f5efd8eb67db6d63a1fcd0990f4f7a21f90 100644 --- a/inmost_options_cmake.h +++ b/inmost_options_cmake.h @@ -15,6 +15,7 @@ #cmakedefine USE_PARTITIONER_PARMETIS #cmakedefine USE_SOLVER +#cmakedefine USE_SOLVER_MONDRIAAN #cmakedefine USE_SOLVER_METIS #cmakedefine USE_SOLVER_PETSC #cmakedefine USE_SOLVER_TRILINOS diff --git a/iterator.cpp b/iterator.cpp index 294839209eebff747741ba3cf51c7de76529dc8f..e3d24b9947fb8b36ad02c598c117b3d115e54a96 100644 --- a/iterator.cpp +++ b/iterator.cpp @@ -50,7 +50,7 @@ namespace INMOST else break; } if( num == 5 && id > 0 ) id = 1; - return ComposeHandle(num,id); + return ComposeHandleNum(num,id); } HandleType Mesh::PrevHandle(HandleType h) const @@ -61,7 +61,7 @@ namespace INMOST { if( id < 0 ) num = 4; - else return ComposeHandle(ElementNum(MESH),0); + else return ComposeHandleNum(ElementNum(MESH),0); } while( num >= 0 ) { @@ -74,7 +74,7 @@ namespace INMOST else break; } if( num < 0 ) return InvalidHandle(); - return ComposeHandle(num,id); + return ComposeHandleNum(num,id); } HandleType Mesh::NextHandle(HandleType h, ElementType etype) const @@ -102,7 +102,7 @@ namespace INMOST else break; } if( num == 5 && id > 0 ) id = 1; - return ComposeHandle(num,id); + return ComposeHandleNum(num,id); } HandleType Mesh::PrevHandle(HandleType h, ElementType etype) const @@ -127,7 +127,7 @@ namespace INMOST } if( stop ) return InvalidHandle(); } - else return ComposeHandle(ElementNum(MESH),0); + else return ComposeHandleNum(ElementNum(MESH),0); } while( num >= 0 ) { @@ -150,7 +150,7 @@ namespace INMOST else break; } if( num < 0 ) return InvalidHandle(); - return ComposeHandle(num,id); + return ComposeHandleNum(num,id); } Storage::integer Mesh::FirstLocalID(ElementType etype) const diff --git a/mesh.cpp b/mesh.cpp index c84cdfea71b8fcc43f46d6456ddd809a52a5875c..82f0a8b97bd700e63734554ae5dc8fe0112009a3 100644 --- a/mesh.cpp +++ b/mesh.cpp @@ -440,7 +440,7 @@ namespace INMOST #endif //USE_MPI #if defined(USE_PARALLEL_WRITE_TIME) FinalizeFile(); - out_time.close(); + for(size_t q = 0; q < allocated_meshes.size(); ++q) if (allocated_meshes[q] == this) allocated_meshes[q] = NULL; @@ -862,7 +862,7 @@ namespace INMOST Node Mesh::CreateNode(const real * coords) { integer id = TieElement(0); - HandleType h = ComposeHandle(0,id); + HandleType h = ComposeHandleNum(0,id); SetGeometricType(h,Element::Vertex); real * v = static_cast(MGetDenseLink(h,CoordsTag())); for(integer i = 0; i < dim; i++) v[i] = coords[i]; @@ -912,7 +912,7 @@ namespace INMOST if (test != InvalidHandle()) return std::make_pair(Edge(this,test),false); } integer id = TieElement(1); - he = ComposeHandle(1,id); + he = ComposeHandleNum(1,id); for(ElementArray::size_type i = 0; i < nodes.size(); i++) { Element::adj_type & hc = HighConn(nodes.at(i)); @@ -1004,7 +1004,7 @@ namespace INMOST if (test != InvalidHandle()) return std::make_pair(Face(this,test),false); } integer id = TieElement(2); - he = ComposeHandle(2,id); + he = ComposeHandleNum(2,id); for(ElementArray::size_type i = 0; i < f_edges.size(); i++) { Element::adj_type & hc = HighConn(f_edges.at(i)); @@ -1374,7 +1374,7 @@ namespace INMOST if (test != InvalidHandle()) return std::make_pair(Cell(this,test),false); } integer id = TieElement(3); - he = ComposeHandle(3,id); + he = ComposeHandleNum(3,id); for(ElementArray::size_type i = 0; i < c_faces.size(); i++) { Element::adj_type & hc = HighConn(c_faces.at(i)); @@ -1476,7 +1476,7 @@ namespace INMOST if( e->GetName() == name ) return std::make_pair(e->self(),false); } - HandleType he = ComposeHandle(4,TieElement(4)); + HandleType he = ComposeHandleNum(4,TieElement(4)); bulk_array set_name = BulkArrayDV(he,SetNameTag()); set_name.resize(static_cast(name.size())); memcpy(set_name.data(),name.c_str(),name.size()); @@ -1498,7 +1498,7 @@ namespace INMOST Storage::integer j = 0; for(Storage::integer k = 0; k < NodeLastLocalID(); ++k) if( isValidElement(0,k) ) { - memcpy(temp.data()+j,MGetDenseLink(ComposeHandle(0,k),CoordsTag()),sizeof(Storage::real)*dims); + memcpy(temp.data()+j,MGetDenseLink(ComposeHandleNum(0,k),CoordsTag()),sizeof(Storage::real)*dims); j+=dims; } @@ -1507,7 +1507,7 @@ namespace INMOST j = 0; for(Storage::integer k = 0; k < NodeLastLocalID(); ++k) if( isValidElement(0,k) ) { - memcpy(MGetDenseLink(ComposeHandle(0,k),CoordsTag()),temp.data()+j,sizeof(Storage::real)*dims); + memcpy(MGetDenseLink(ComposeHandleNum(0,k),CoordsTag()),temp.data()+j,sizeof(Storage::real)*dims); j+=dims; } dim = dims; @@ -1536,7 +1536,7 @@ namespace INMOST back_links[etypenum][ADDR] = -1; empty_space[etypenum].push_back(ADDR); empty_links[etypenum].push_back(ID); - //REPORT_VAL("destroyed",ComposeHandle(etypenum,ID) << " " << etypenum << " " << ADDR << " " << ID); + //REPORT_VAL("destroyed",ComposeHandleNum(etypenum,ID) << " " << etypenum << " " << ADDR << " " << ID); } } @@ -1576,7 +1576,7 @@ namespace INMOST new_size = GetArrayCapacity(etypenum); if( new_size != old_size ) ReallocateData(etypenum,new_size); back_links[etypenum][ADDR] = ID; - last_created = ComposeHandle(etypenum,ID); + last_created = ComposeHandleNum(etypenum,ID); //REPORT_VAL("created",last_created << " " << etypenum << " " << ADDR << " " << ID); } return ID; diff --git a/mesh_file.cpp b/mesh_file.cpp index 8a13c84786d4bf24827e31d51c92719956b722fb..96e11f7e62199a56fadb3636ec077ec3e19496c7 100644 --- a/mesh_file.cpp +++ b/mesh_file.cpp @@ -2946,6 +2946,7 @@ read_elem_num_link: } else if(LFile.find(".pmf") != std::string::npos) //this is inner parallel/platform mesh format { + //std::cout << "parallel strategy " << parallel_strategy << " file strategy " << parallel_file_strategy << std::endl; io_converter iconv; io_converter uconv; REPORT_STR("start load pmf"); @@ -2958,13 +2959,16 @@ read_elem_num_link: #if defined(USE_MPI_FILE) if( parallel_file_strategy == 1 ) { + REPORT_STR("strategy 1"); int ierr; - std::string buffer; + std::vector buffer; MPI_File fh; MPI_Status stat; - ierr = MPI_File_open(GetCommunicator(),const_cast(File.c_str()),MPI_MODE_RDONLY,MPI_INFO_NULL,&fh); - if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),-1); + REPORT_MPI(ierr = MPI_File_open(GetCommunicator(),const_cast(File.c_str()),MPI_MODE_RDONLY,MPI_INFO_NULL,&fh)); + if( ierr != MPI_SUCCESS ) REPORT_MPI(MPI_Abort(GetCommunicator(),-1)); INMOST_DATA_ENUM_TYPE numprocs = GetProcessorsNumber(), recvsize = 0, mpirank = GetProcessorRank(); + REPORT_VAL("number of processors",numprocs); + REPORT_VAL("rank of processor", mpirank); std::vector recvsizes; if( mpirank == 0 ) //the alternative is to read alltogether { @@ -2973,32 +2977,44 @@ read_elem_num_link: buffer.resize(3); //ierr = MPI_File_read_all(fh,&buffer[0],3,MPI_CHAR,&stat); - ierr = MPI_File_read_shared(fh,&buffer[0],3,MPI_CHAR,&stat); - if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),-1); + REPORT_MPI(ierr = MPI_File_read_shared(fh,&buffer[0],3,MPI_CHAR,&stat)); + if( ierr != MPI_SUCCESS ) REPORT_MPI(MPI_Abort(GetCommunicator(),-1)); if( static_cast(buffer[0]) != INMOST::INMOSTFile ) throw BadFile; - header << buffer.c_str()+1; + header.write(&buffer[1],2); uconv.read_iByteOrder(header); uconv.read_iByteSize(header); + REPORT_VAL("integer_byte_order",uconv.str_iByteOrder(uconv.get_iByteOrder())); + REPORT_VAL("integer_byte_size",(int)uconv.get_iByteSize()); + + buffer.resize(uconv.get_source_iByteSize()); //ierr = MPI_File_read_all(fh,&buffer[0],buffer.size(),MPI_CHAR,&stat); - ierr = MPI_File_read_shared(fh,&buffer[0],static_cast(buffer.size()),MPI_CHAR,&stat); + REPORT_MPI(ierr = MPI_File_read_shared(fh,&buffer[0],static_cast(buffer.size()),MPI_CHAR,&stat)); if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),-1); - header << buffer; + header.write(&buffer[0],buffer.size()); uconv.read_iValue(header,datanum); + REPORT_VAL("number of data entries",datanum); + buffer.resize(datanum*uconv.get_source_iByteSize()); std::vector datasizes(datanum); //ierr = MPI_File_read_all(fh,&buffer[0],buffer.size(),MPI_CHAR,&stat); - ierr = MPI_File_read_shared(fh,&buffer[0],static_cast(buffer.size()),MPI_CHAR,&stat); - if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),-1); - - header << buffer; - for(k = 0; k < datanum; k++) uconv.read_iValue(header,datasizes[k]); + REPORT_MPI(ierr = MPI_File_read_shared(fh,&buffer[0],static_cast(buffer.size()),MPI_CHAR,&stat)); + if( ierr != MPI_SUCCESS ) REPORT_MPI(MPI_Abort(GetCommunicator(),-1)); + INMOST_DATA_ENUM_TYPE datasum = 0; + header.write(&buffer[0],buffer.size()); + for(k = 0; k < datanum; k++) + { + uconv.read_iValue(header,datasizes[k]); + REPORT_VAL("size of data entry " << k,datasizes[k]); + datasum += datasizes[k]; + } + REPORT_VAL("total size",datasum); // use this commented code when all processors read the file alltogether through MPI_File_read_all //{ @@ -3031,92 +3047,134 @@ read_elem_num_link: recvsizes.resize(numprocs,0); if( datanum <= recvsizes.size() ) { + REPORT_STR("number of processors is greater or equal then number of data entries"); for(k = 0; k < datanum; k++) recvsizes[k] = datasizes[k]; } else { + REPORT_STR("number of processors is less then number of data entries - accumulating data"); chunk = static_cast(floor(static_cast(datanum)/static_cast(recvsizes.size()))); + REPORT_VAL("chunk size" ,chunk); pos = 0; for(k = 0; k < recvsizes.size()-1; k++) + { for(q = 0; q < chunk; q++) recvsizes[k] += datasizes[pos++]; + + REPORT_VAL("recv on " << k, recvsizes[k]); + } for(k = pos; k < datanum; k++) recvsizes[recvsizes.size()-1] += datasizes[k]; + REPORT_VAL("recv on " << recvsizes.size()-1, recvsizes[recvsizes.size()-1]); } } else recvsizes.resize(1,0); //protect from dereferencing null - ierr = MPI_Scatter(&recvsizes[0],1,INMOST_MPI_DATA_ENUM_TYPE,&recvsize,1,INMOST_MPI_DATA_ENUM_TYPE,0,GetCommunicator()); - if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),-1); + REPORT_MPI(ierr = MPI_Scatter(&recvsizes[0],1,INMOST_MPI_DATA_ENUM_TYPE,&recvsize,1,INMOST_MPI_DATA_ENUM_TYPE,0,GetCommunicator())); + if( ierr != MPI_SUCCESS ) REPORT_MPI(MPI_Abort(GetCommunicator(),-1)); + REPORT_VAL("read on current processor",recvsize); buffer.resize(std::max(1u,recvsize)); //protect from dereferencing null { - ierr = MPI_File_read_ordered(fh,&buffer[0],static_cast(recvsize),MPI_CHAR,&stat); - if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),-1); + REPORT_MPI(ierr = MPI_File_read_ordered(fh,&buffer[0],static_cast(recvsize),MPI_CHAR,&stat)); + if( ierr != MPI_SUCCESS ) REPORT_MPI(MPI_Abort(GetCommunicator(),-1)); in.write(&buffer[0],recvsize); } - ierr = MPI_File_close(&fh); - if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),-1); + REPORT_MPI(ierr = MPI_File_close(&fh)); + if( ierr != MPI_SUCCESS ) REPORT_MPI(MPI_Abort(GetCommunicator(),-1)); } else #endif { + + REPORT_STR("strategy 0"); int ierr; - std::string buffer, local_buffer; + std::vector buffer, local_buffer; INMOST_DATA_ENUM_TYPE recvsize; - INMOST_DATA_ENUM_TYPE numprocs = GetProcessorsNumber(); + INMOST_DATA_ENUM_TYPE numprocs = GetProcessorsNumber(),mpirank = GetProcessorRank(); std::vector recvsizes(numprocs,0); std::vector sendcnts(numprocs), displs(numprocs); - if( GetProcessorRank() == 0 ) //zero reads everything + REPORT_VAL("number of processors",numprocs); + REPORT_VAL("rank of processor", mpirank); + if( mpirank == 0 ) //zero reads everything { std::fstream fin(File.c_str(),std::ios::in | std::ios::binary); fin.get(token); if( token != INMOST::INMOSTFile ) throw BadFile; uconv.read_iByteOrder(fin); uconv.read_iByteSize(fin); + + REPORT_VAL("file position",fin.tellg()); + + REPORT_VAL("integer_byte_order",uconv.str_iByteOrder(uconv.get_iByteOrder())); + REPORT_VAL("integer_byte_size",(int)uconv.get_iByteSize()); + + INMOST_DATA_ENUM_TYPE datanum,k,q,datasum = 0,chunk,pos; uconv.read_iValue(fin,datanum); + + REPORT_VAL("number of data entries",datanum); + REPORT_VAL("file position",fin.tellg()); + std::vector datasizes(datanum); for(k = 0; k < datanum; k++) { uconv.read_iValue(fin,datasizes[k]); + REPORT_VAL("size of data entry " << k,datasizes[k]); datasum += datasizes[k]; } + + REPORT_VAL("file position",fin.tellg()); + { buffer.resize(datasum); fin.read(&buffer[0],buffer.size()); + } + REPORT_VAL("file position",fin.tellg()); + REPORT_VAL("total size",datasum); fin.close(); if( datanum <= recvsizes.size() ) { + REPORT_STR("number of processors is greater or equal then number of data entries"); for(k = 0; k < datanum; k++) recvsizes[k] = datasizes[k]; } else { + REPORT_STR("number of processors is less then number of data entries - accumulating data"); chunk = static_cast(floor(static_cast(datanum)/static_cast(recvsizes.size()))); + REPORT_VAL("chunk size" ,chunk); pos = 0; for(k = 0; k < recvsizes.size()-1; k++) + { for(q = 0; q < chunk; q++) recvsizes[k] += datasizes[pos++]; + REPORT_VAL("recv on " << k, recvsizes[k]); + } for(k = pos; k < datanum; k++) recvsizes[recvsizes.size()-1] += datasizes[k]; + REPORT_VAL("recv on " << recvsizes.size()-1, recvsizes[recvsizes.size()-1]); } displs[0] = 0; sendcnts[0] = static_cast(recvsizes[0]); + REPORT_VAL("disp on "<<0,displs[0]); + REPORT_VAL("send on "<<0,sendcnts[0]); for(k = 1; k < numprocs; k++) { sendcnts[k] = static_cast(recvsizes[k]); displs[k] = sendcnts[k-1]+displs[k-1]; + REPORT_VAL("disp on "< datasizes(datanum); for(k = 0; k < datanum; k++) { uconv.read_iValue(fin,datasizes[k]); + REPORT_VAL("size of data entry " << k,datasizes[k]); datasum += datasizes[k]; } + + REPORT_VAL("file position",fin.tellg()); + + REPORT_VAL("total size",datasum); { - std::string buffer; + std::vector buffer; buffer.resize(datasum); fin.read(&buffer[0],buffer.size()); in.write(&buffer[0],buffer.size()); + REPORT_VAL("output position",in.tellg()); } + + REPORT_VAL("file position",fin.tellg()); + fin.close(); } //std::fstream in(File.c_str(),std::ios::in | std::ios::binary); std::vector tags; + std::vector tags_defined; + std::vector tags_sparse; std::vector old_nodes; std::vector new_nodes; std::vector new_edges; @@ -3179,6 +3264,7 @@ read_elem_num_link: while (in >> token) { + REPORT_VAL("output position, loop",in.tellg()); if( !start ) { if( token != INMOST::INMOSTFile ) throw BadFile; //check that this is valid file @@ -3187,6 +3273,8 @@ read_elem_num_link: REPORT_STR("File chunk start read"); //~ std::cout << "start read" << std::endl; tags.clear(); + tags_sparse.clear(); + tags_defined.clear(); old_nodes.clear(); old_nodes.resize(NumberOfNodes()); { @@ -3272,6 +3360,8 @@ read_elem_num_link: new_sets.clear(); new_sets.resize(header[5]); tags.resize(header[6]); + tags_sparse.resize(header[6]); + tags_defined.resize(header[6]); //~ if( static_cast(header[7]) == Mesh::Parallel && m_state != Mesh::Parallel) //~ SetCommunicator(INMOST_MPI_COMM_WORLD); myprocs.push_back(header[8]); @@ -3292,13 +3382,25 @@ read_elem_num_link: in.read(name, namesize); assert(namesize < 4096); name[namesize] = '\0'; + REPORT_VAL("tag name",name); in.get(datatype); + REPORT_VAL("tag data type",DataTypeName(static_cast(datatype))); in.get(sparsemask); in.get(definedmask); + //for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype) ) + //{ + // if( etype & definedmask ) REPORT_VAL("defined on",ElementTypeName(etype)); + // if( etype & sparsemask ) REPORT_VAL("sparse on",ElementTypeName(etype)); + //} uconv.read_iValue(in,datalength); + REPORT_VAL("length",datalength); tags[i] = CreateTag(std::string(name),static_cast(datatype), static_cast(definedmask), static_cast(sparsemask),datalength); + tags_defined[i] = static_cast(definedmask); + tags_sparse[i] = static_cast(sparsemask); + + REPORT_VAL("output position, tag " << i,in.tellg()); } } else if (token == INMOST::NodeHeader) @@ -3325,6 +3427,8 @@ read_elem_num_link: if( find == -1 ) new_nodes[i] = CreateNode(coords)->GetHandle(); else new_nodes[i] = old_nodes[find]; } + + REPORT_VAL("output position, nodes",in.tellg()); } else if (token == INMOST::EdgeHeader) { @@ -3345,6 +3449,8 @@ read_elem_num_link: new_edges[i] = CreateEdge(sub_elements).first->GetHandle(); sub_elements.clear(); } + + REPORT_VAL("output position, edges",in.tellg()); } else if (token == INMOST::FaceHeader) { @@ -3365,6 +3471,8 @@ read_elem_num_link: new_faces[i] = CreateFace(sub_elements).first->GetHandle(); sub_elements.clear(); } + + REPORT_VAL("output position, faces",in.tellg()); } else if (token == INMOST::CellHeader) { @@ -3393,6 +3501,8 @@ read_elem_num_link: sub_elements.clear(); suggest_nodes.clear(); } + + REPORT_VAL("output position, cells",in.tellg()); } else if (token == INMOST::ESetHeader) { @@ -3481,6 +3591,8 @@ read_elem_num_link: lc[j] = new_sets[GetHandleID(lc[j])]; } } + + REPORT_VAL("output position, sets",in.tellg()); } else if (token == INMOST::MeshDataHeader) { @@ -3511,20 +3623,28 @@ read_elem_num_link: Tag * jt = &tags[j]; for(ElementType etype = NODE; etype <= MESH; etype = etype << 1) { - if( jt->isDefined(etype) ) + if( etype & tags_defined[j] ) { + REPORT_VAL("defined on",ElementTypeName(etype)); INMOST_DATA_ENUM_TYPE q, cycle_end, etypenum = ElementNum(etype); cycle_end = elem_sizes[etypenum]; - bool sparse = jt->isSparse(etype); + REPORT_VAL("cycle end",cycle_end); + bool sparse = false; + if( etype & tags_sparse[j] ) sparse = true; INMOST_DATA_ENUM_TYPE tagsize = jt->GetSize(), recsize = tagsize, lid; INMOST_DATA_ENUM_TYPE k; DataType data_type = jt->GetDataType(); if( sparse ) { + REPORT_VAL("sparse on",ElementTypeName(etype)); uconv.read_iValue(in,q); if( in.eof() ) std::cout << __FILE__ << ":" << __LINE__ << " Unexpected end of file! " << tags[j].GetTagName() << " " << ElementTypeName(etype) << " " << (sparse? "sparse" : "dense") << std::endl; } else q = 0; + + REPORT_VAL("data type",DataTypeName(data_type)); + REPORT_VAL("tag size",tagsize); + while(q != cycle_end) { HandleType he = elem_links[etypenum][q]; @@ -3586,7 +3706,11 @@ read_elem_num_link: } } } + + REPORT_VAL("output position, tag data " << j,in.tellg()); } + + REPORT_VAL("output position, tag data",in.tellg()); if( verbosity > 0 ) std::cout << "Finished reading data" << std::endl; REPORT_STR("EndOfData"); } diff --git a/mesh_parallel.cpp b/mesh_parallel.cpp index b0710cf2e17cc5ff6c58ff545a28bd4cdddab6fb..b7f7ad260fc39ff7f1e3e558901cf5d881fefa51 100644 --- a/mesh_parallel.cpp +++ b/mesh_parallel.cpp @@ -76,7 +76,7 @@ namespace INMOST { if( size ) { - const Storage::integer * recv = static_cast(static_cast(data)); + const Storage::integer * recv = static_cast(static_cast(data)); Storage::integer_array arr = e->IntegerArray(tag); arr.push_back(recv[0]); } @@ -87,20 +87,26 @@ namespace INMOST { if( size ) { - bool flag = true; - const Storage::integer * recv = static_cast(static_cast(data)); - Storage::integer_array arr = e->IntegerArray(tag); - for(Storage::integer_array::iterator it = arr.begin(); it != arr.end(); it+=2) - if( *it == recv[0] ) - { - flag = false; - break; - } - if( flag ) - { - arr.push_back(recv[0]); - arr.push_back(recv[1]); - } + for(INMOST_DATA_ENUM_TYPE k = 0; k < size; k+=2 ) + { + bool flag = true; + const Storage::integer * recv = static_cast(static_cast(data))+k; + Storage::integer_array arr = e->IntegerArray(tag); + + for(Storage::integer_array::iterator it = arr.begin(); it != arr.end(); it+=2) + { + if( *it == recv[0] ) + { + flag = false; + break; + } + } + if( flag ) + { + arr.push_back(recv[0]); + arr.push_back(recv[1]); + } + } } } @@ -302,9 +308,9 @@ namespace INMOST #endif//USE_MPI shift += start; for(const HandleType * it = set; it != set+n; ++it) - if( GetStatus(*it) != Element::Owned && (define_sparse || HaveData(*it,num_tag))) Integer(*it,num_tag) = shift++; + if( GetStatus(*it) == Element::Owned && (define_sparse || HaveData(*it,num_tag))) Integer(*it,num_tag) = shift++; for(const HandleType * it = set; it != set+n; ++it) - if( GetStatus(*it) != Element::Shared && (define_sparse || HaveData(*it,num_tag))) Integer(*it,num_tag) = shift++; + if( GetStatus(*it) == Element::Shared && (define_sparse || HaveData(*it,num_tag))) Integer(*it,num_tag) = shift++; ExchangeData(num_tag,mask,0); ret = shift; #if defined(USE_MPI) @@ -326,10 +332,10 @@ namespace INMOST #endif//USE_MPI shift += start; for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++) - if( it->GetStatus() != Element::Owned && (define_sparse || it->HaveData(num_tag)) ) + if( it->GetStatus() == Element::Owned && (define_sparse || it->HaveData(num_tag)) ) it->Integer(num_tag) = shift++; for(Mesh::iteratorElement it = BeginElement(mask); it != EndElement(); it++) - if( it->GetStatus() != Element::Shared && (define_sparse || it->HaveData(num_tag)) ) + if( it->GetStatus() == Element::Shared && (define_sparse || it->HaveData(num_tag)) ) it->Integer(num_tag) = shift++; ExchangeData(num_tag,mask,0); ret = shift; @@ -625,7 +631,7 @@ namespace INMOST #if defined(USE_PARALLEL_WRITE_TIME) void Mesh::Enter() { tab++; } void Mesh::Exit() {tab--; } - std::fstream & Mesh::WriteTab(std::fstream & f) + std::ostream & Mesh::WriteTab(std::ostream & f) { for(int i = 0; i < tab; i++) f << " "; @@ -637,15 +643,20 @@ namespace INMOST } void Mesh::FinalizeFile() { - if( tab > 1 ) + //std::stringstream str; + if( tab > 1 ) { - REPORT_STR("Finalizing file with nonempty function stack - probably due to error."); + out_time << "\n";// << std::endl; } - while(tab > 1) + while(tab > 1) { - EXIT_FUNC_DIE(); + out_time << "\n\n";// << std::endl; + Exit(); } out_time << "" << std::endl; + //out_time << str; + //out_time.flush(); + out_time.close(); } #endif //USE_PARALLEL_WRITE_TIME @@ -1005,7 +1016,7 @@ namespace INMOST //~ #else std::vector sendsizeall(mpisize*2); int pack_size2 = 0; - unsigned long usend[2] = {static_cast(sendsize),static_cast(pack_real.size())}; + unsigned usend[2] = {static_cast(sendsize),static_cast(pack_real.size())}; MPI_Pack_size(2,MPI_UNSIGNED,comm,&pack_size2); for(dynarray::size_type k = 0; k < procs.size(); k++) { @@ -1160,6 +1171,8 @@ namespace INMOST REPORT_VAL("type",ElementTypeName(current_mask)); + //int owned_elems = 0; + //int shared_elems = 0; int owner; Element::Status estat; @@ -1173,10 +1186,19 @@ namespace INMOST { p.clear(); p.push_back(mpirank); + //++owned_elems; } - else p.replace(p.begin(),p.end(),result.begin(),result.end()); + else + { + p.replace(p.begin(),p.end(),result.begin(),result.end()); + //if( result.size() == 1 && result[0] == mpirank ) + // ++owned_elems; + //else ++shared_elems; + } } time = Timer() - time; + //REPORT_VAL("predicted owned elements",owned_elems); + //REPORT_VAL("predicted shared elements",shared_elems); REPORT_STR("Predict processors for elements"); REPORT_VAL("time",time); @@ -1184,6 +1206,7 @@ namespace INMOST time = Timer(); //Initialize mapping that helps get local id by global id std::vector > mapping; + REPORT_VAL("mapping type",ElementTypeName(current_mask >> 1)); for(Mesh::iteratorElement it = BeginElement(current_mask >> 1); it != EndElement(); it++) { mapping.push_back(std::make_pair(it->GlobalID(),it->LocalID())); @@ -1191,6 +1214,7 @@ namespace INMOST if( !mapping.empty() ) std::sort(mapping.begin(),mapping.end(),MappingComparator()); time = Timer() - time; + REPORT_VAL("mapping size",mapping.size()) REPORT_STR("Compute global to local indexes mapping"); REPORT_VAL("time",time); //Initialize arrays @@ -1227,6 +1251,8 @@ namespace INMOST elements[m].push_back(*it); } } + REPORT_VAL("for processor",*p); + REPORT_VAL("gathered elements",elements[m].size()); message_send[0] = static_cast(message_send.size()); MPI_Pack_size(static_cast(message_send.size()),MPI_INT,comm,&sendsize); send_buffs[m].first = *p; @@ -1264,6 +1290,7 @@ namespace INMOST } if( pos == -1 ) throw Impossible; MPI_Unpack(&recv_buffs[*qt].second[0],static_cast(recv_buffs[*qt].second.size()),&position,&size,1,MPI_INT,comm); + REPORT_VAL("unpacked message size",size-1); message_recv[pos].resize(size-1); MPI_Unpack(&recv_buffs[*qt].second[0],static_cast(recv_buffs[*qt].second.size()),&position,&message_recv[pos][0],static_cast(message_recv[pos].size()),MPI_INT,comm); } @@ -1301,7 +1328,7 @@ namespace INMOST break; } int find_local_id = mapping[find].second; - sub_elements.push_back(ComposeHandle(current_mask >> 1, find_local_id)); + sub_elements.push_back(ComposeHandle(PrevElementType(current_mask), find_local_id)); } if( flag ) @@ -1311,14 +1338,32 @@ namespace INMOST remote_elements.push_back(e); } } + REPORT_VAL("number of unpacked remote elements",remote_elements.size()); + if( !remote_elements.empty() ) + { + REPORT_VAL("first",remote_elements.front()); + REPORT_VAL("first type",ElementTypeName(GetHandleElementType(remote_elements.front()))); + REPORT_VAL("last",remote_elements.back()); + REPORT_VAL("last type",ElementTypeName(GetHandleElementType(remote_elements.back()))); + } std::sort(remote_elements.begin(),remote_elements.end()); + REPORT_VAL("original elements size",elements[m].size()); + if( !elements[m].empty() ) + { + REPORT_VAL("first",elements[m].front()); + REPORT_VAL("first type",ElementTypeName(GetHandleElementType(elements[m].front()))); + REPORT_VAL("last",elements[m].back()); + REPORT_VAL("last type",ElementTypeName(GetHandleElementType(elements[m].back()))); + } std::sort(elements[m].begin(),elements[m].end()); element_set result; element_set::iterator set_end; result.resize(elements[m].size()); set_end = std::set_difference(elements[m].begin(),elements[m].end(),remote_elements.begin(),remote_elements.end(), result.begin()); - result.resize(set_end-result.begin()); + result.resize(set_end-result.begin()); + + REPORT_VAL("set difference size",result.size()); //elements in result are wrongly marked as ghost for(element_set::iterator qt = result.begin(); qt != result.end(); qt++) @@ -1799,6 +1844,9 @@ namespace INMOST Tag tag_skin = CreateTag("TEMPORARY_COMPUTE_SHARED_SKIN_SET",DATA_INTEGER,FACE,FACE); REPORT_STR("filling neighbouring cell's global identificators for faces") +#if defined(USE_PARALLEL_WRITE_TIME) + std::map numfacesperproc; +#endif for(iteratorFace it = BeginFace(); it != EndFace(); it++) { @@ -1807,30 +1855,68 @@ namespace INMOST { Storage::integer_array arr = it->IntegerArray(tag_skin); ElementArray adj = it->getAdjElements(CELL); + //REPORT_STR("face " << it->LocalID() << " global " << it->GlobalID() << " type " << Element::StatusName(estat)); for(ElementArray::iterator jt = adj.begin(); jt != adj.end(); jt++) { arr.push_back(jt->GlobalID()); //identificator of the cell arr.push_back(jt->IntegerDF(tag_owner)); //owner of the cell + //REPORT_STR("id " << jt->GlobalID() << " owner " << jt->IntegerDF(tag_owner)); } } +#if defined(USE_PARALLEL_WRITE_TIME) + ++numfacesperproc[it->IntegerDF(tag_owner)]; +#endif } + REPORT_STR("number of ghosted or shared faces"); +#if defined(USE_PARALLEL_WRITE_TIME) + for(std::map::iterator it = numfacesperproc.begin(); it != numfacesperproc.end(); ++it) + { + REPORT_VAL("processor",it->first); + REPORT_VAL("skin faces",it->second); + } +#endif - REPORT_STR("exchanging cell's global identificators information") + REPORT_STR("reducing cell's global identificators information") ReduceData(tag_skin,FACE,0,UnpackSkin); + + + REPORT_STR("exchanging cell's global identificators information") + +#if defined(USE_PARALLEL_WRITE_TIME) + //for(iteratorFace it = BeginFace(); it != EndFace(); it++) + //{ + // Element::Status estat1 = GetStatus(*it), estat2; + // if( estat1 == Element::Owned ) continue; + // Storage::integer_array skin_data = it->IntegerArray(tag_skin); + // REPORT_STR("face " << it->LocalID() << " global " << it->GlobalID() << " type " << Element::StatusName(estat1)); + // for(Storage::integer_array::iterator kt = skin_data.begin(); kt != skin_data.end(); kt+=2) + // { + // REPORT_STR("id " << *kt << " owner " << *(kt+1)); + // } + //} +#endif + ExchangeData(tag_skin,FACE,0); - REPORT_STR("exchanging cell's global identificators information") + + REPORT_STR("synchornization done") proc_elements skin_faces; #if defined(USE_PARALLEL_WRITE_TIME) - std::map numfacesperproc; + numfacesperproc.clear(); #endif for(iteratorFace it = BeginFace(); it != EndFace(); it++) { bool flag = false; Element::Status estat1 = GetStatus(*it), estat2; if( estat1 == Element::Owned ) continue; + //Storage::integer_array skin_data = it->IntegerArray(tag_skin); + //REPORT_STR("face " << it->LocalID() << " global " << it->GlobalID() << " type " << Element::StatusName(estat1)); + //for(Storage::integer_array::iterator kt = skin_data.begin(); kt != skin_data.end(); kt+=2) + //{ + // REPORT_STR("id " << *kt << " owner " << *(kt+1)); + //} Cell c1 = it->BackCell(); Cell c2 = it->FrontCell(); if( !c1.isValid() && !c2.isValid() ) continue; //no cells per face - skip hanging face @@ -1856,6 +1942,7 @@ namespace INMOST if( flag ) { + //printf("hello!\n"); for(int i = 0; i < 2; i++) //assert checks that there are two cells { Storage::integer owner = skin_data[i*2+1]; //cell owner @@ -1963,7 +2050,7 @@ namespace INMOST for(element_set::iterator it = all_visited.begin(); it != all_visited.end(); it++) { Storage::integer_array os = IntegerArray(*it,on_skin); - for(Storage::integer_array::iterator p = os.begin(); p != os.end(); p++) + for(Storage::integer_array::iterator p = os.begin(); p != os.end(); p++) { if( *p != mpirank ) { @@ -2019,12 +2106,18 @@ namespace INMOST for(eit = elements[i].begin(); eit != elements[i].end(); eit++) if( (!select || GetMarker(*eit,select)) && HaveData(*eit,tag) ) { + //REPORT_STR("element type " << ElementTypeName(GetHandleElementType(*eit)) << " global id " << Integer(*eit,GlobalIDTag())); array_size_send.push_back(static_cast(eit-elements[i].begin())); array_size_send[count]++; INMOST_DATA_ENUM_TYPE s = GetDataSize(*eit,tag); INMOST_DATA_ENUM_TYPE had_s = static_cast(array_data_send.size()); array_data_send.resize(had_s+s*tag.GetBytesSize()); GetData(*eit,tag,0,s,&array_data_send[had_s]); + //REPORT_VAL("size",s); + //for(int qq = 0; qq < s; ++qq) + //{ + // REPORT_VAL("value " << qq, (*(Storage::integer *)&array_data_send[had_s+qq*tag.GetBytesSize()])); + //} if( size == ENUMUNDEF ) array_size_send.push_back(s); ++total_packed; } @@ -2033,11 +2126,12 @@ namespace INMOST { for(eit = elements[i].begin(); eit != elements[i].end(); eit++) if( !select || GetMarker(*eit,select) ) { + // REPORT_STR("element type " << ElementTypeName(*eit) << " global id " << Integer(*eit,GlobalIDTag())); INMOST_DATA_ENUM_TYPE s = GetDataSize(*eit,tag); INMOST_DATA_ENUM_TYPE had_s = static_cast(array_data_send.size()); array_data_send.resize(had_s+s*tag.GetBytesSize()); GetData(*eit,tag,0,s,&array_data_send[had_s]); - if( size == ENUMUNDEF ) array_size_send.push_back(s); + if( size == ENUMUNDEF ) array_size_send.push_back(s); ++total_packed; } } @@ -2112,7 +2206,13 @@ namespace INMOST for(unsigned j = 0; j < count; j++) { eit = elements[i].begin() + array_size_recv[k++]; + //REPORT_STR("element type " << ElementTypeName(GetHandleElementType(*eit)) << " global id " << Integer(*eit,GlobalIDTag())); assert( !select || GetMarker(*eit,select) ); //if fires then very likely that marker was not synchronized + //REPORT_VAL("size",array_size_recv[k]); + //for(int qq = 0; qq < array_size_recv[k]; ++qq) + //{ + // REPORT_VAL("value " << qq, (*(Storage::integer *)&array_data_recv[pos+qq*tag.GetBytesSize()])); + //} op(tag,Element(this,*eit),&array_data_recv[pos],array_size_recv[k]); pos += array_size_recv[k]*tag.GetBytesSize(); ++k; @@ -3980,7 +4080,7 @@ namespace INMOST void Mesh::ExchangeGhost(Storage::integer layers, ElementType bridge) { - //printf("called exchange ghost with %d bridge %s\n",layers,ElementTypeName(bridge)); + //printf("%d called exchange ghost with layers %d bridge %s\n",GetProcessorRank(), layers,ElementTypeName(bridge)); if( m_state == Serial ) return; ENTER_FUNC(); #if defined(USE_MPI) @@ -4007,6 +4107,7 @@ namespace INMOST { proc_elements shared_skin = ComputeSharedSkinSet(bridge); + //printf("%d shared skin size %d\n",GetProcessorRank(),shared_skin.size()); time = Timer(); for(Storage::integer_array::iterator p = procs.begin(); p != procs.end(); p++) { diff --git a/partitioner.cpp b/partitioner.cpp index ddd081cd44a2e2a62263588481d09f972a465274..67e950d0ab5f482d6b2b38d698e1a6a45d45d7db 100644 --- a/partitioner.cpp +++ b/partitioner.cpp @@ -346,9 +346,17 @@ namespace INMOST { queue.push_back(*it); m->Integer(queue.back(),index) = visited++; + + if( visited % chunk == 0 ) + { + queue.clear(); + break; + } + } } - if( visited % chunk == 0 ) queue.clear(); + //this may not fire properly + //if( visited % chunk == 0 ) queue.clear(); } } m->DeleteTag(order); diff --git a/solver.cpp b/solver.cpp index 0999f122fe55569512bbb9158fc719cf78d3934e..0b04c0210f521ac7e2106dce65d899e8781bb680 100644 --- a/solver.cpp +++ b/solver.cpp @@ -1,3 +1,4 @@ +#define _CRT_SECURE_NO_WARNINGS #include "inmost_solver.h" #if defined(USE_SOLVER) #include "solver_petsc.h" @@ -1337,37 +1338,47 @@ namespace INMOST (void)argv; if( database != NULL ) { - std::fstream file(database,std::ios::in); - char str[4096]; - while( !file.eof() && file.getline(str,4096) ) - { - int k = 0, l; - for(k = 0; k < (int)strlen(str); ++k) - { - if( str[k] == ':' ) break; - } - if( k == strlen(str) ) continue; //invalid line - for(l = 0; l < k; ++l) str[l] = tolower(str[l]); - l = k+1; - while(l < (int)strlen(str) && isspace(str[l]) ) ++l; - if( l == strlen(str) ) continue; //skip empty entry - if( !strncmp(str,"petsc",k) ) - petsc_database_file = std::string(str+l); - else if( !strncmp(str,"trilinos_ifpack",k) ) - trilinos_ifpack_database_file = std::string(str+l); - else if( !strncmp(str,"trilinos_aztec",k) ) - trilinos_aztec_database_file = std::string(str+l); - else if( !strncmp(str,"trilinos_ml",k) ) - trilinos_ml_database_file = std::string(str+l); - else if( !strncmp(str,"trilinos_belos",k) ) - trilinos_belos_database_file = std::string(str+l); - else if( !strncmp(str,"ani",k) ) - ani_database_file = std::string(str+l); - else if( !strncmp(str,"fcbiilu2",k) ) - fcbiilu2_database_file = std::string(str+l); - else if( !strncmp(str,"k3biilu2",k) ) - k3biilu2_database_file = std::string(str+l); - } + FILE * f = fopen(database,"r"); + if( f ) + { + //std::fstream file(database,std::ios::in); + char str[4096]; + //while( !file.eof() && file.getline(str,4096) ) + while( !feof(f) && fgets(str,4096,f) ) + { + int k = 0, l; + for(k = 0; k < (int)strlen(str); ++k) + { + if( str[k] == ':' ) break; + } + if( k == strlen(str) ) continue; //invalid line + for(l = 0; l < k; ++l) str[l] = tolower(str[l]); + l = (int)strlen(str)-1; // Right-trim string + while(l > 0 && isspace(str[l]) ) --l; + str[l+1] = 0; + l = k+1; + while(l < (int)strlen(str) && isspace(str[l]) ) ++l; + if( l == strlen(str) ) continue; //skip empty entry + if( !strncmp(str,"petsc",k) ) + petsc_database_file = std::string(str+l); + else if( !strncmp(str,"trilinos_ifpack",k) ) + trilinos_ifpack_database_file = std::string(str+l); + else if( !strncmp(str,"trilinos_aztec",k) ) + trilinos_aztec_database_file = std::string(str+l); + else if( !strncmp(str,"trilinos_ml",k) ) + trilinos_ml_database_file = std::string(str+l); + else if( !strncmp(str,"trilinos_belos",k) ) + trilinos_belos_database_file = std::string(str+l); + else if( !strncmp(str,"ani",k) ) + ani_database_file = std::string(str+l); + else if( !strncmp(str,"fcbiilu2",k) ) + fcbiilu2_database_file = std::string(str+l); + else if( !strncmp(str,"k3biilu2",k) ) + k3biilu2_database_file = std::string(str+l); + } + //file.close(); + fclose(f); + } } //std::cout << "PETSc \"" << petsc_database_file << "\"" << std::endl; //std::cout << "Trilinos_Ifpack \"" << trilinos_ifpack_database_file << "\"" << std::endl; diff --git a/solver_bcgsl.hpp b/solver_bcgsl.hpp index 8d4da97e898d993470aa767f2ce874c2eb2b760f..3ff9981ea3cb7f1182429053f548271410417e7c 100644 --- a/solver_bcgsl.hpp +++ b/solver_bcgsl.hpp @@ -46,13 +46,15 @@ namespace INMOST shell u(pu,n*n); shell v(pv,n*n); shell w(pw,n); - std::copy(a.begin(),a.end(),u.begin()); + //std::copy(a.begin(),a.end(),u.begin()); //memcpy(u,a,sizeof(INMOST_DATA_REAL_TYPE)*n*n); int flag, i, its, j, jj, k, l, nm; INMOST_DATA_REAL_TYPE c, f, h, s, x, y, z; INMOST_DATA_REAL_TYPE anorm = 0.0, g = 0.0, scale = 0.0; - static std::vector rv1(n); + dynarray rv1; + rv1.resize(n); // Householder reduction to bidiagonal form + for (i = 0; i < n*n; ++i) u[i] = a[i]; for (i = 0; i < n; i++) { // left-hand reduction @@ -624,13 +626,13 @@ namespace INMOST } #endif INMOST_DATA_ENUM_TYPE i = 0; - bool halt = false; if( last_resid < atol || last_resid < rtol*resid0 ) { reason = "initial solution satisfy tolerances"; - goto exit; + halt = true; } + #if defined(USE_OMP) #pragma omp parallel #endif @@ -826,17 +828,24 @@ namespace INMOST sigma[j-1] = sigma_sum; } #if defined(CONVEX_COMBINATION) - INMOST_DATA_REAL_TYPE lagrangian = 0.0; - for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++) lagrangian += tau[j+size*j]; - sigma[l] = lagrangian; - tau[(l+1)*(l+1)-1] = 0.0; - for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++) +#if defined(USE_OMP) +#pragma omp single +#endif { - tau[l + j*(l+1)] = -lagrangian; - tau[l*(l+1)+j] = lagrangian; + INMOST_DATA_REAL_TYPE lagrangian = 0.0; + for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++) lagrangian += tau[j+size*j]; + sigma[l] = lagrangian; + tau[(l+1)*(l+1)-1] = 0.0; + for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++) + { + tau[l + j*(l+1)] = -lagrangian; + tau[l*(l+1)+j] = lagrangian; + } } #endif - info->Integrate(tau,size*size+size); //sigma is updated with tau + info->Integrate(tau,(l+2)*(l+1)); //sigma is updated with tau + //info->Integrate(tau,size*size); + //info->Integrate(sigma,size); #if defined(USE_OMP) #pragma omp single @@ -1136,7 +1145,6 @@ namespace INMOST i++; } } -exit: if (prec != NULL) { prec->Solve(SOL, r_tilde); //undo right preconditioner @@ -1277,6 +1285,10 @@ exit: if (prec != NULL)prec->ReplaceRHS(RHS); info->GetLocalRegion(info->GetRank(),vlocbeg,vlocend); info->GetVectorRegion(vbeg,vend); + ivbeg = vbeg; + ivend = vend; + ivlocbeg = vlocbeg; + ivlocend = vlocend; std::copy(RHS.Begin(),RHS.End(),r.Begin()); { Alink->MatVec(-1,SOL,1,r); //global multiplication, r probably needs an update @@ -1296,10 +1308,6 @@ exit: } last_resid = resid = resid0 = sqrt(resid); last_it = 0; - ivbeg = vbeg; - ivend = vend; - ivlocbeg = vlocbeg; - ivlocend = vlocend; #if defined(REPORT_RESIDUAL) if( info->GetRank() == 0 ) { diff --git a/tests/solver_test000/CMakeLists.txt b/tests/solver_test000/CMakeLists.txt index c5af2675133fe9d599011fddd8d61a74afc69de1..57dcef533ad3a6c538a25a9ced7554cbec7befd9 100644 --- a/tests/solver_test000/CMakeLists.txt +++ b/tests/solver_test000/CMakeLists.txt @@ -32,6 +32,10 @@ if(USE_SOLVER_METIS) message("linking solver_test000 with Metis") target_link_libraries(solver_test000 ${METIS_LIBRARIES}) endif() +if(USE_SOLVER_MONDRIAAN) + message("linking solver_test000 with Mondriaan") + target_link_libraries(solver_test000 ${MONDRIAAN_LIBRARIES}) +endif() add_test(NAME solver_test000_serial_inner_ilu2 COMMAND $ 0 0) diff --git a/tests/solver_test001/CMakeLists.txt b/tests/solver_test001/CMakeLists.txt index 6a602a26a242261c939d393f9575748dec3ac4c8..903092c5440ac71e1bf4e19b272c352a0fb2f27f 100644 --- a/tests/solver_test001/CMakeLists.txt +++ b/tests/solver_test001/CMakeLists.txt @@ -30,6 +30,10 @@ if(USE_SOLVER_METIS) message("linking solver_test001 with Metis") target_link_libraries(solver_test001 ${METIS_LIBRARIES}) endif() +if(USE_SOLVER_MONDRIAAN) + message("linking solver_test001 with Mondriaan") + target_link_libraries(solver_test001 ${MONDRIAAN_LIBRARIES}) +endif() set(TOKAMAK_TESTS utm300 diff --git a/tests/solver_test002/CMakeLists.txt b/tests/solver_test002/CMakeLists.txt index e519cf90fb673a9a5903e8a3c18ef09ec539e269..7ad5d2551e39f22d286b3ef9bfef56fa5b92f268 100644 --- a/tests/solver_test002/CMakeLists.txt +++ b/tests/solver_test002/CMakeLists.txt @@ -33,6 +33,10 @@ if(USE_SOLVER) message("linking solver_test002 with Metis") target_link_libraries(solver_test002 ${METIS_LIBRARIES}) endif() + if(USE_SOLVER_MONDRIAAN) + message("linking solver_test002 with Mondriaan") + target_link_libraries(solver_test002 ${MONDRIAAN_LIBRARIES}) + endif() endif() add_test(NAME solver_test002_serial_inner_ilu2 COMMAND $ 0 20)