Commit 07f25538 by Kirill Terekhov

Fix possible bug, new grid tool

Possible fix in packtagdata/unpacktagdata - mismatch between amount of
calls to MPI_Pack and MPI_Unpack

GridTool for Kmeans clustering, should eventually move to Partitioner
class
parent 58fab942
......@@ -16,6 +16,7 @@ add_executable(Move move.cpp)
add_executable(Convert convert.cpp)
add_executable(SameMeshDifference difference_same.cpp)
add_executable(MeshDifference difference_map.cpp)
add_executable(Kmeans kmeans.cpp)
add_library(FractureLib fracture.cpp fracture.h)
target_link_libraries(FixFaults inmost)
......@@ -185,6 +186,16 @@ if(USE_MPI)
endif(USE_MPI)
install(TARGETS MeshDifference EXPORT inmost-targets RUNTIME DESTINATION bin)
target_link_libraries(Kmeans inmost)
if(USE_MPI)
message("linking Kmeans with MPI")
target_link_libraries(Kmeans ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(Kmeans PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Kmeans EXPORT inmost-targets RUNTIME DESTINATION bin)
set_property(TARGET FractureLib PROPERTY PUBLIC_HEADER "${CMAKE_CURRENT_SOURCE_DIR}/fracture.h")
......
#include "inmost.h"
using namespace INMOST;
//todo: want to separate all faces into those that share edges
//see todo: inside text
double dot(const double a[3], const double b[3])
{
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
struct coords
{
double v[3];
double operator [] (int i) const {return v[i];}
double & operator [] (int i) {return v[i];}
int compare(const double bv[3]) const
{
const double eps = 1.0e-11;//pow(10,-14);
for (int i = 0; i < 3; i++)
{
if (fabs(v[i] - bv[i]) > eps)
{
if (v[i] > bv[i])
return 1;
else
return -1;
}
}
return 0;
}
bool operator <(const coords & b) const { return compare(b.v) < 0; }
coords() {v[0] = v[1] = v[2] = 0;}
coords(double x, double y, double z)
{
v[0] = x;
v[1] = y;
v[2] = z;
}
coords(double _v[3]) { memcpy(v, _v, sizeof(double)* 3); }
coords(const coords &b) { memcpy(v, b.v, sizeof(double)* 3); }
coords & operator =(coords const & b) { memmove(v, b.v, sizeof(double)* 3); return *this; }
coords operator -(const coords & b) const { return coords(v[0] - b.v[0], v[1] - b.v[1], v[2] - b.v[2]); }
coords operator +(const coords & b) const { return coords(v[0] + b.v[0], v[1] + b.v[1], v[2] + b.v[2]); }
coords operator /(double l) const {return coords(v[0]/l, v[1]/l, v[2]/l);}
};
double len2(const coords & p)
{
return p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
}
struct Point
{
coords c;
int cluster;
int position;
int node;
void set_cluster(int c, int p)
{
cluster = c;
position = p;
}
coords & get_center()
{
return c;
}
};
struct Cluster
{
std::vector<int> points;
std::vector<int> emptyp;
coords center;
int add_position(int p)
{
if( emptyp.empty() )
{
points.push_back(p);
return points.size()-1;
}
else
{
int pos = emptyp.back();
points[pos] = p;
emptyp.pop_back();
return pos;
}
}
void rem_position(int k)
{
points[k] = -1;
emptyp.push_back(k);
}
void set_center(coords c)
{
center = c;
}
int num_points()
{
return (int)(points.size() - emptyp.size());
}
coords & get_center() {return center;}
};
int main(int argc, char ** argv)
{
if (argc < 2)
{
std::cout << "Usage: " << argv[0] << " mesh clusters [clusters=10] [max_iterations=10] [mesh_out=grid.vtk]" << std::endl;
return -1;
}
int K = 10;
if( argc > 2 ) K = atoi(argv[2]);
int max_iterations = 10;
if( argc > 3 ) max_iterations = atoi(argv[3]);
if( K == 0 )
std::cout << "Problem with number of clusters argument " << argv[2] << std::endl;
if( max_iterations == 0 )
std::cout << "Problem with max iterations argument " << argv[3] << std::endl;
std::string grid_out = "grid.vtk";
if (argc > 4) grid_out = std::string(argv[4]);
Mesh m;
m.Load(argv[1]);
std::cout << "Start:" << std::endl;
std::cout << "Cells: " << m.NumberOfCells() << std::endl;
std::cout << "Faces: " << m.NumberOfFaces() << std::endl;
std::cout << "Edges: " << m.NumberOfEdges() << std::endl;
double tt = Timer();
//There seems to be a problem with the mesh
int fixed = 0;
for (Mesh::iteratorFace f = m.BeginFace(); f != m.EndFace(); ++f)
{
if (f->FixNormalOrientation())
fixed++;
}
std::cout << "Time to fix normals: " << Timer() - tt << std::endl;
std::cout << "Total face normals fixed: " << fixed << std::endl;
int total_points = m.NumberOfCells();
std::vector< Cluster > clusters(K);
std::vector< Point > points(total_points);
int k = 0;
double cnt[3];
for(Mesh::iteratorCell n = m.BeginCell(); n != m.EndCell(); ++n)
{
n->Centroid(cnt);
points[k].c = coords(cnt);
points[k].cluster = -1;
points[k].position = -1;
points[k].node = n->LocalID();
k++;
}
// choose K distinct values for the centers of the clusters
{
std::vector<int> prohibited_indexes;
for(int i = 0; i < K; i++)
{
while(true)
{
int index_point = rand() % total_points;
if(find(prohibited_indexes.begin(), prohibited_indexes.end(),
index_point) == prohibited_indexes.end())
{
prohibited_indexes.push_back(index_point);
int pos = clusters[i].add_position(index_point);
points[index_point].set_cluster(i,pos);
clusters[i].set_center(points[index_point].c);
break;
}
}
}
}
int iter = 1;
while(true)
{
bool done = true;
// associates each point to the nearest center
for(int i = 0; i < total_points; i++)
{
int id_old_cluster = points[i].cluster;
int id_nearest_center = -1;
double lmin = 1.0e+100;
for(int j = 0; j < K; ++j)
{
double l = len2(points[i].get_center() - clusters[j].get_center());
if( l < lmin )
{
lmin = l;
id_nearest_center = j;
}
}
if(id_old_cluster != id_nearest_center)
{
if(id_old_cluster != -1)
clusters[id_old_cluster].rem_position(points[i].position);
int pos = clusters[id_nearest_center].add_position(i);
points[i].set_cluster(id_nearest_center,pos);
done = false;
}
}
// recalculating the center of each cluster
for(int i = 0; i < K; i++)
{
int total_points_cluster = clusters[i].num_points();
coords sum(0,0,0);
if(total_points_cluster > 0)
{
for(int p = 0; p < (int)clusters[i].points.size(); p++)
if( clusters[i].points[p] != -1 )
sum = sum + points[clusters[i].points[p]].get_center();
clusters[i].set_center(sum / (double) total_points_cluster);
}
}
if(done == true || iter >= max_iterations)
{
std::cout << "Break in iteration " << iter << "\n\n";
break;
}
else std::cout << "Iteration " << iter << std::endl;
iter++;
}
// shows elements of clusters
TagInteger mat = m.CreateTag("MATERIAL",DATA_INTEGER,CELL,NONE,1);
for(int i = 0; i < K; i++)
{
for(int j = 0; j < (int)clusters[i].points.size(); ++j)
if( clusters[i].points[j] != -1 )
mat[m.CellByLocalID(points[clusters[i].points[j]].node)] = i+1;
}
m.Save(grid_out);
return 0;
}
......@@ -2296,7 +2296,7 @@ namespace INMOST
ElementType pack_types[2] = {NONE,NONE};
element_set::const_iterator eit;
buffer_type array_data_send;
std::vector<INMOST_DATA_ENUM_TYPE> array_size_send(2);
std::vector<INMOST_DATA_ENUM_TYPE> array_size_send;
array_data_send.reserve(4096);
array_size_send.reserve(4096);
unsigned int size = tag.GetSize();
......@@ -2381,8 +2381,9 @@ namespace INMOST
}
REPORT_VAL("total packed records",total_packed);
}
array_size_send[0] = static_cast<INMOST_DATA_ENUM_TYPE>(array_size_send.size()-2);
array_size_send[1] = static_cast<INMOST_DATA_ENUM_TYPE>(array_data_send.size());
INMOST_DATA_ENUM_TYPE size_send, data_send;
size_send = static_cast<INMOST_DATA_ENUM_TYPE>(array_size_send.size());
data_send = static_cast<INMOST_DATA_ENUM_TYPE>(array_data_send.size());
REPORT_VAL("tag defined on",static_cast<int>(pack_types[0]));
REPORT_VAL("tag sparse on",static_cast<int>(pack_types[1]));
REPORT_VAL("size_size",array_size_send[0]);
......@@ -2390,10 +2391,14 @@ namespace INMOST
int buffer_size = 0,position = static_cast<int>(buffer.size()),temp,bytes;
bytes = tag.GetPackedBytesSize();
MPI_Pack_size(2 ,INMOST_MPI_DATA_BULK_TYPE,comm,&temp); buffer_size+= temp;
MPI_Pack_size(1 ,INMOST_MPI_DATA_ENUM_TYPE,comm,&temp); buffer_size+= temp;
MPI_Pack_size(1 ,INMOST_MPI_DATA_ENUM_TYPE,comm,&temp); buffer_size+= temp;
MPI_Pack_size(static_cast<INMOST_MPI_SIZE>(array_size_send.size()) ,INMOST_MPI_DATA_ENUM_TYPE,comm,&temp); buffer_size+= temp;
MPI_Pack_size(static_cast<INMOST_MPI_SIZE>(array_data_send.size()/bytes),tag.GetBulkDataType() ,comm,&temp); buffer_size+= temp;
buffer.resize(position+buffer_size);
MPI_Pack(pack_types,2,INMOST_MPI_DATA_BULK_TYPE,&buffer[0],static_cast<INMOST_MPI_SIZE>(buffer.size()),&position,comm);
MPI_Pack(&size_send,1,INMOST_MPI_DATA_ENUM_TYPE,&buffer[0],static_cast<INMOST_MPI_SIZE>(buffer.size()),&position,comm);
MPI_Pack(&data_send,1,INMOST_MPI_DATA_ENUM_TYPE,&buffer[0],static_cast<INMOST_MPI_SIZE>(buffer.size()),&position,comm);
if( !array_size_send.empty() ) MPI_Pack(&array_size_send[0],static_cast<INMOST_MPI_SIZE>(array_size_send.size()),INMOST_MPI_DATA_ENUM_TYPE,&buffer[0],static_cast<INMOST_MPI_SIZE>(buffer.size()),&position,comm);
if( !array_data_send.empty() ) MPI_Pack(&array_data_send[0],static_cast<INMOST_MPI_SIZE>(array_data_send.size()/bytes),tag.GetBulkDataType(),&buffer[0],static_cast<INMOST_MPI_SIZE>(buffer.size()),&position,comm);
buffer.resize(position);
......@@ -2468,7 +2473,7 @@ namespace INMOST
array_size_recv.resize(size_recv);
array_data_recv.resize(data_recv);
if( !array_size_recv.empty() ) MPI_Unpack(&buffer[0],static_cast<INMOST_MPI_SIZE>(buffer.size()),&position,&array_size_recv[0],static_cast<INMOST_MPI_SIZE>(array_size_recv.size()),INMOST_MPI_DATA_ENUM_TYPE,comm);
if( !array_data_recv.empty() )
if( !array_data_recv.empty() )
{
int bytes = tag.GetPackedBytesSize();
REPORT_VAL("occupied by type",bytes);
......
......@@ -24,7 +24,7 @@ using namespace INMOST;
#define REORDER_RCM
//#define REORDER_NNZ
#if defined(USE_SOLVER_METIS)
//#define REORDER_METIS_ND
#define REORDER_METIS_ND
#endif
#if defined(USE_SOLVER_MONDRIAAN)
//#define REORDER_MONDRIAAN
......@@ -37,13 +37,13 @@ static bool rescale_b = true;
static bool allow_pivot = true;
#define ESTIMATOR
//#define ESTIMATOR_REFINE
#define ESTIMATOR_REFINE
//#define PREMATURE_DROPPING
//#define EQUALIZE_1NORM
#define EQUALIZE_2NORM
//#define EQUALIZE_IDOMINANCE
//#define EQUALIZE_2NORM
#define EQUALIZE_IDOMINANCE
#define PIVOT_THRESHOLD
#define PIVOT_THRESHOLD_VALUE 1.0e-12
......@@ -52,7 +52,7 @@ static bool allow_pivot = true;
#define DIAGONAL_PERTURBATION_ABS 1.0e-10
#define ILUC2
#define ILUC2_SCHUR
//#define TRACK_DIAGONAL
#define TRACK_DIAGONAL
//#define PIVOT_COND_DEFAULT 0.1/tau
#define PIVOT_COND_DEFAULT 1.0e+2
#define PIVOT_DIAG_DEFAULT 1.0e+5
......
......@@ -26,7 +26,7 @@ using namespace INMOST;
#define REORDER_RCM
//#define REORDER_NNZ
#if defined(USE_SOLVER_METIS)
//#define REORDER_METIS_ND
#define REORDER_METIS_ND
#endif
#if defined(USE_SOLVER_MONDRIAAN)
//#define REORDER_MONDRIAAN
......@@ -45,15 +45,15 @@ using namespace INMOST;
//#define EQUALIZE_IDOMINANCE
#define PIVOT_THRESHOLD
#define PIVOT_THRESHOLD_VALUE 1.0e-6
#define DIAGONAL_PERTURBATION
#define PIVOT_THRESHOLD_VALUE 1.0e-9
//#define DIAGONAL_PERTURBATION
#define DIAGONAL_PERTURBATION_REL 1.0e-7
#define DIAGONAL_PERTURBATION_ABS 1.0e-10
//#define DIAGONAL_PIVOT //probably there is some bug
#define DIAGONAL_PIVOT_TAU 0.01
//#define DIAGONAL_PIVOT_COND 100
#define ILUC2
#define TRACK_DIAGONAL
//#define TRACK_DIAGONAL
#if defined(DIAGONAL_PIVOT) && defined(DIAGONAL_PIVOT_TAU) && !defined(TRACK_DIAGONAL)
#define TRACK_DIAGONAL
......@@ -443,6 +443,7 @@ using namespace INMOST;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> Ubeg(mobeg, moend,EOL), Lbeg(mobeg, moend,EOL), Bbeg(mobeg,moend,EOL);
INMOST_DATA_REAL_TYPE NuU = 1, NuL = 1, NuD, NuU_max = 1.0, NuL_max = 1.0;
INMOST_DATA_REAL_TYPE NuUsqrt = 1, NuLsqrt = 1;
#if defined(ESTIMATOR)
//supplimentary data structures for condition estimates of L^{-1}, U^{-1}
INMOST_DATA_REAL_TYPE mup, mum, smup, smum, NuL1 = 1, NuL2 = 1, NuU1 = 1, NuU2 = 1;
......@@ -2754,6 +2755,9 @@ swap_algorithm:
#endif
NuL_max = std::max(NuL,NuL_max);
NuU_max = std::max(NuU,NuU_max);
NuUsqrt = sqrt(NuU);
NuLsqrt = sqrt(NuL);
}
#endif //ESTIMATOR
......@@ -2794,7 +2798,7 @@ swap_algorithm:
{
u = fabs(LineValuesL[Li]);
//if( !(u == u) ) std::cout << __FILE__ << ":" << __LINE__ << " nan " << std::endl;
if (u*NuL /*NuU_old*/ > tau) //apply dropping
if (u*NuL /*NuU_old*/ > tau) //apply dropping
LU_Entries.push_back(Sparse::Row::make_entry(Li, LineValuesL[Li]));
#if defined(ILUC2)
else if (u*NuL /*NuU_old*/ > tau2)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment