Commit 9a92cf91 authored by Igor Konshin's avatar Igor Konshin

Merge branch 'master' of https://github.com/INMOST-DEV/INMOST

Conflicts:
	solver.cpp
parents b2d17c6f 8cfc4bfa
......@@ -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)
......
......@@ -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.
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.
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
......@@ -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()
......
......@@ -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()
......
......@@ -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
......
......@@ -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()
......
......@@ -538,6 +538,7 @@ class incident_matrix
}
return success;
}
/*
Storage::real compute_measure(dynarray<T,64> & 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<T,64>::iterator it = min_loop.begin(); it != min_loop.end(); ++it)
for(typename dynarray<T,64>::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();
......
......@@ -11,6 +11,17 @@
#include <algorithm>
#include <stdarg.h>
#include "my_glut.h"
#include <iomanip>
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<Element> 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<float> ticks; //ticks from 0 to 1 for each color
std::vector<color_t> 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<float>::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<double> verts;
std::vector<color_t> 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<face2gl> & set)
{
static std::vector<face2gl> 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<float>(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<face2gl> all_boundary;
std::vector<face2gl> clip_boundary;
void draw_faces_nc(std::vector<face2gl> & 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<face2gl> & set, int highlight = -1)
{
glBegin(GL_TRIANGLES);
......@@ -190,6 +486,13 @@ void draw_faces(std::vector<face2gl> & set, int highlight = -1)
}
}
void draw_faces_alpha(std::vector<face2gl> & 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<face2gl> & set, int highlight = -1)
{
glBegin(GL_LINES);
......@@ -205,6 +508,14 @@ void draw_edges(std::vector<face2gl> & set, int highlight = -1)
}
void draw_faces_interactive_nc(std::vector<face2gl> & 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<face2gl> & set)
{
glBegin(GL_TRIANGLES);
......@@ -212,6 +523,13 @@ void draw_faces_interactive(std::vector<face2gl> & set)
glEnd();
}
void draw_faces_interactive_alpha(std::vector<face2gl> & 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<face2gl> & 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<point_t> points;
Mesh * m;
void radix_sort_dist(std::vector<point_t> & set)
{
static std::vector<point_t> 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<Face> faces = it->getFaces();