Commit b1919dfd authored by SilverLife's avatar SilverLife

Fix depth bug in Octree example

parent f87c239c
......@@ -20,6 +20,7 @@ int draw_faces = 1;
int draw_ghost = 1;
int draw_in_motion = 0; // redraw mode. true - always, false - by space press.
int draw_sem = 0;
int all_cells_count;
int cur_cell = 0; bool draw_all_cells = true;
int cur_face = 0; bool draw_all_faces = true;
......@@ -29,10 +30,10 @@ int from_file = 0;
int redist_after_amr = 0;
// Variables for refine and coarse
int refine_depth = 1;
int refine_depth = 2;
double mx = 0, my = 0;
double rmx = 0, rmy = 0;
double base_radius = 0.02;
double base_radius = 0.01;
int action = 0;
int log_level = 0;
......@@ -91,6 +92,14 @@ void dump_to_vtk()
cout << "Process " << ::rank << ": dumped mesh to file" << endl;
}
void print_all_cells(grid* g)
{
for(Mesh::iteratorCell it = thegrid.mesh->BeginCell(); it != thegrid.mesh->EndCell(); it++)
{
cout << ::rank << "-" << it->Integer(g->c_tags.base_id) << endl;
}
}
/// Function provided to octgrid algorithm. Defines transformation from grid to grid for draw.
void transformation(double xyz[3])
{
......@@ -231,8 +240,10 @@ void mpi_draw()
else if (it->GetStatus() == Element::Shared) glColor3f(1.0, 1.0 ,0.0);
else glColor3f(0.0, 0.0 ,1.0);
if (!draw_all_cells && cur_cell % thegrid.mesh->NumberOfCells() != c++) continue;
if (!draw_all_cells && cur_cell % thegrid.mesh->NumberOfCells() != c++) continue;
// Draw faces
ElementArray<Face> faces = it->getFaces();
for (ElementArray<Face>::iterator f = faces.begin(); f != faces.end(); f++)
......@@ -277,7 +288,6 @@ void mpi_draw()
for (int j = 1; j < drawing_faces_n[i]; j++)
{
if (!draw_all_faces && cur_face % drawing_faces_n[i] != cf++) continue;
if (!draw_all_faces) cout << drawing_faces[i][j].nodes_count << endl;
glBegin(GL_POLYGON);
for (int k = 0; k < drawing_faces[i][j].nodes_count; k++)
{
......@@ -519,11 +529,13 @@ void prepare_to_correct_brothers()
/// Redistribute grid by partitioner
void redistribute(int type)
{
thegrid.mesh->RemoveGhost();
LOG(2,"Process " << ::rank << ": redistribute. Cells: " << thegrid.mesh->NumberOfCells())
Partitioner * part = new Partitioner(thegrid.mesh);
// Specify the partitioner
type = 0;
type = 1;
if (type == 0) part->SetMethod(Partitioner::Parmetis, Partitioner::Partition);
if (type == 1) part->SetMethod(Partitioner::Parmetis, Partitioner::Repartition);
if (type == 2) part->SetMethod(Partitioner::Parmetis, Partitioner::Refine);
......@@ -541,6 +553,7 @@ void redistribute(int type)
}
delete part;
thegrid.mesh->RemoveGhost();
correct_brothers(&thegrid,size,::rank, 2);
try
......@@ -872,11 +885,73 @@ void NotMainProcess()
}
}
void parse_arguments(int argc, char** argv, int* n, double* R)
{
if (argc < 2) return;
string str;
string str1, str2;
for (int i = 1; i < argc; i++)
{
str = argv[i];
size_t pos = str.find('=');
if (pos == string::npos)
{
cout << "Invalid argument: " << str << endl;
continue;
}
str1 = str.substr(0, pos);
str2 = str.substr(pos + 1);
if (str1 == "-n")
{
for (int j = 0; j < 3; j++)
{
pos = str2.find('x');
if (j < 2 && pos == string::npos)
{
cout << "Invalid command: " << str2 << endl;
break;
}
str1 = str2.substr(0, pos);
n[j] = atoi(str1.c_str());
str2 = str2.substr(pos + 1);
}
}
else if (str1 == "-r")
{
*R = atof(str2.c_str());
}
else
{
cout << "Invalid command: " << str1 << endl;
}
}
}
void print_help()
{
cout << "Example of Octree refine on redistributed grid" << endl;
cout << "Command arguments:" << endl;
cout << " -n=20x20x1 - grid size" << endl;
cout << " -R=0.02 - refine radius" << endl;
cout << endl;
cout << "Hotkeys:" << endl;
cout << " Space - refine grid around mouse cursor" << endl;
cout << " r - redraw grid" << endl;
cout << " f - dump grid to file (see grids folder)" << endl;
cout << " x - redistribute grid" << endl;
}
int main(int argc, char ** argv)
{
int i;
int n[3] = {10,10,1};
thegrid.transformation = transformation;
thegrid.rev_transformation = rev_transformation;
thegrid.cell_should_split = cell_should_split;
......@@ -884,6 +959,10 @@ int main(int argc, char ** argv)
Mesh::Initialize(&argc,&argv);
MPI_Comm_rank(MPI_COMM_WORLD, &::rank);
if (::rank == 0) print_help();
parse_arguments(argc, argv, n, &base_radius);
all_cells_count = n[0]*n[1]*n[2] * 2;
gridInit(&thegrid,n);
Partitioner::Initialize(&argc,&argv);
......
......@@ -2,6 +2,7 @@
#include <math.h>
#include <iostream>
#include <iomanip>
#include <map>
using namespace std;
double epsilon = 0.0001;
......@@ -16,6 +17,8 @@ int default_cell_should_split(struct grid * g, int cell) { (void) g; (void) cell
void mark_cell(struct grid* g, Cell cell);
void cellCoarse(struct grid* g, Cell cell);
void resolve_connections(grid* g);
inline bool equal(double a, double b);
int global = 0;
void dump_to_vtk(grid* g)
{
......@@ -47,8 +50,30 @@ void init_mesh(struct grid* g)
g->c_tags.chld_num = g->mesh->CreateTag("chld_num",DATA_INTEGER,CELL,false,1);
g->c_tags.par_chld_nums = g->mesh->CreateTag("par_chld_nums",DATA_INTEGER,CELL,false,8);
g->c_tags.to_split = g->mesh->CreateTag("to_split",DATA_INTEGER,CELL,false,1);
g->c_tags.base_id = g->mesh->CreateTag("base_id" ,DATA_INTEGER,CELL,false,1);
}
int check_all_edges(grid* g, int n)
{
int rank = g->mesh->GetProcessorRank(); // Get the rank of the current process
int was = 0;
for(Mesh::iteratorEdge edge = g->mesh->BeginEdge(); edge != g->mesh->EndEdge(); edge++)
{
Node n1 = edge->getBeg();
Node n2 = edge->getEnd();
int diff = 0;
if (!equal(n1.RealArray(g->mesh->CoordsTag())[0], n2.RealArray(g->mesh->CoordsTag())[0])) diff++;
if (!equal(n1.RealArray(g->mesh->CoordsTag())[1], n2.RealArray(g->mesh->CoordsTag())[1])) diff++;
if (!equal(n1.RealArray(g->mesh->CoordsTag())[2], n2.RealArray(g->mesh->CoordsTag())[2])) diff++;
if (diff > 1)
{
print_edge(g,edge->getAsEdge());
was = 1;
}
}
return was;
}
/// Debug
void print_node_center(struct grid* g, Node node)
{
......@@ -252,6 +277,7 @@ void gridInit(struct grid * g, int n[3])
g->mesh = new Mesh();
g->mesh->SetDimensions(3);
g->last_base_id = 0;
g->mesh->SetCommunicator(INMOST_MPI_COMM_WORLD);
MPI_Comm_set_errhandler(INMOST_MPI_COMM_WORLD,MPI_ERRORS_RETURN);
......@@ -360,6 +386,7 @@ void gridInit(struct grid * g, int n[3])
c.Integer(g->c_tags.floor) = k;
c.Integer(g->c_tags.level) = 0;
c.Integer(g->c_tags.chld_num) = -1;
c.Integer(g->c_tags.base_id) = g->last_base_id++ + (rank * n[0]*n[1]*n[2]);;
for (int i = 0; i < 8; i++)
c.IntegerArrayDF(g->c_tags.par_chld_nums)[i] = -1;
......@@ -476,6 +503,7 @@ Cell recreate_cell(struct grid* g, Cell cell, ElementArray<Face> faces)
res.Integer(g->c_tags.level) = cell.Integer(g->c_tags.level);
res.Integer(g->c_tags.to_split) = cell.Integer(g->c_tags.to_split);
res.Integer(g->c_tags.base_id) = cell.Integer(g->c_tags.base_id);
g->mesh->Destroy(cell);
return res;
......@@ -587,6 +615,27 @@ double dot(double* a, double* b, int n = 3)
return res;
}
void order_nodes_in_face_b(struct grid* g, ElementArray<Node>* nodes)
{
/*
int rank = g->mesh->GetProcessorRank();
for (int i = 0; i < nodes->size() - 1; i++)
{
double* nc = (*nodes)[i].RealArray(g->mesh->CoordsTag()).
double min_dist = dist(nc,(*nodes)[i+1].RealArray(g->mesh->CoordsTag()));
int min_j = i+1;
for (int j = i+2; j < nodes->size(); j++)
{
double* ac = (*nodes)[j].RealArray(g->mesh->CoordsTag());
if (dist(nc,ac) < min_dist)
{
min_dist = dist(nc,ac);
min_j = j;
}
}
}
*/
}
/// Reorder nodes with this order
/// 1 2
/// 0 3
......@@ -604,7 +653,7 @@ void order_nodes_in_face_a(struct grid* g, ElementArray<Node>* nodes)
if (i < nodes->size() - 1)
{
for (int j = 0;j < 3;j++)
if ((*nodes)[i].RealArray(g->mesh->CoordsTag())[j] != (*nodes)[i+1].RealArray(g->mesh->CoordsTag())[j]) xyz[j] = 0;
if (!equal((*nodes)[i].RealArray(g->mesh->CoordsTag())[j],(*nodes)[i+1].RealArray(g->mesh->CoordsTag())[j])) xyz[j] = 0;
}
cen[0] += (*nodes)[i].RealArray(g->mesh->CoordsTag())[0];
cen[1] += (*nodes)[i].RealArray(g->mesh->CoordsTag())[1];
......@@ -613,7 +662,7 @@ void order_nodes_in_face_a(struct grid* g, ElementArray<Node>* nodes)
cen[0] /= nodes->size();
cen[1] /= nodes->size();
cen[2] /= nodes->size();
if (xyz[0] == 0 && xyz[1] == 0) { I[0] = 0; I[1] = 1; }
else if (xyz[1] == 0 && xyz[2] == 0) { I[0] = 1; I[1] = 2; }
else if (xyz[0] == 0 && xyz[2] == 0) { I[0] = 0; I[1] = 2; }
......@@ -706,6 +755,7 @@ ElementArray<Cell> cellSplit_2(struct grid* g, Cell cell)
{
int rank = g->mesh->GetProcessorRank();
int ghost = IS_GHOST(cell);
ElementArray<Edge> edges = cell.getEdges();
for (ElementArray<Edge>::iterator edge = edges.begin(); edge != edges.end(); edge++)
{
......@@ -826,6 +876,7 @@ ElementArray<Cell> cellSplit_2(struct grid* g, Cell cell)
double ndz = dz/2.0;
int level = cell.Integer(g->c_tags.level);
int floor = cell.Integer(g->c_tags.floor);
int base_id = cell.Integer(g->c_tags.base_id);
ElementArray<Node> local_nodes;
ElementArray<Cell> result;
......@@ -849,6 +900,7 @@ ElementArray<Cell> cellSplit_2(struct grid* g, Cell cell)
c.Integer(g->c_tags.level) = level + 1;
c.Integer(g->c_tags.floor) = floor;
c.Integer(g->c_tags.chld_num) = i;
c.Integer(g->c_tags.base_id) = base_id;
bool ready = false;
for (int k = 0; k < 8; k++)
{
......@@ -878,6 +930,7 @@ int get_adj_count(Node node)
/// Unite array children to one cell
Cell cell_unite(struct grid* g, ElementArray<Cell> children, Node center)
{
int rank = g->mesh->GetProcessorRank();
double dx = children[0].RealArrayDF(g->c_tags.side)[0];
double dy = children[0].RealArrayDF(g->c_tags.side)[1];
double dz = children[0].RealArrayDF(g->c_tags.side)[2];
......@@ -886,6 +939,7 @@ Cell cell_unite(struct grid* g, ElementArray<Cell> children, Node center)
double ndz = dz*2.0;
int level = children[0].Integer(g->c_tags.level);
int floor = children[0].Integer(g->c_tags.floor);
int base_id = children[0].Integer(g->c_tags.base_id);
double xyz[3];
xyz[0] = center->RealArray(g->mesh->CoordsTag())[0];
xyz[1] = center->RealArray(g->mesh->CoordsTag())[1];
......@@ -915,7 +969,7 @@ Cell cell_unite(struct grid* g, ElementArray<Cell> children, Node center)
for (ElementArray<Cell>::iterator p = children.begin(); p != children.end(); p++)
{
ElementArray<Node> nodes = p->getNodes();
if (nodes.size() != 8) TSNH
if (nodes.size() != 8) { cout << rank << "|" << nodes.size(); TSNH }
reorder_1(g,&nodes);
int j = 0;
while (j < 5 && verts[p->Integer(g->c_tags.chld_num)][j] != -1)
......@@ -956,7 +1010,7 @@ Cell cell_unite(struct grid* g, ElementArray<Cell> children, Node center)
Face new_face = g->mesh->CreateFace(face_nodes).first;
ElementArray<Face> faces = def_nodes[ x_nodes[i] ].getFaces();
//if (faces.size() != 4) TSNH
if (faces.size() != 4) { cout << faces.size(); TSNH }
if (faces[0].getCells().size() != 1) { cout << faces[0].getCells().size(); TSNH }
if (def_nodes[x_nodes[i]].getCells().size() != 1) TSNH
......@@ -972,7 +1026,7 @@ Cell cell_unite(struct grid* g, ElementArray<Cell> children, Node center)
q->getAsFace() != faces[3]) new_neight_faces.push_back(q->getAsFace());
}
new_neight_faces.push_back(new_face);
if (neight_faces.size() - new_neight_faces.size() != 3) TSNH
if (neight_faces.size() - new_neight_faces.size() != 3) { TSNH }
recreate_cell(g, neight_cell, new_neight_faces);
new_faces.push_back(new_face);
......@@ -1018,6 +1072,7 @@ Cell cell_unite(struct grid* g, ElementArray<Cell> children, Node center)
c.Integer(g->c_tags.floor) = floor;
c.Integer(g->c_tags.level) = level - 1;
c.Integer(g->c_tags.base_id) = base_id;
for (int i = 0; i < 8; i++)
{
......@@ -1165,7 +1220,7 @@ void cellCoarse(struct grid* g, Cell cell)
ElementArray<Cell> children;
ElementArray<Element> adjs = center_node.getAdjElements(CELL);
if (adjs.size() != 8) TSNH
if (adjs.size() != 8) { cout << adjs.size(); TSNH }
for (ElementArray<Element>::iterator adj = adjs.begin(); adj != adjs.end(); adj++)
{
......@@ -1325,13 +1380,76 @@ void print_redist_tag(struct grid* g, int rank)
}
}
void correct_brothers(struct grid* g, int size, int rank, int type)
Face get_bottom_or_upper_face(grid* g, Cell c, bool bottom)
{
ElementArray<Node> nodes = c.getNodes();
if (nodes.size() == 0) TSNH
double minmax_z = nodes[0].RealArray(g->mesh->CoordsTag())[2];
for (int i = 1; i < nodes.size(); i++)
{
if (bottom && nodes[i].RealArray(g->mesh->CoordsTag())[2] < minmax_z
|| !bottom && nodes[i].RealArray(g->mesh->CoordsTag())[2] > minmax_z)
{
minmax_z = nodes[i].RealArray(g->mesh->CoordsTag())[2];
}
}
ElementArray<Face> faces = c.getFaces();
if (faces.size() < 1) TSNH
Face result_face = faces[0];
for (int i = 0; i < faces.size(); i++)
{
ElementArray<Node> f_nodes = faces[i].getNodes();
bool done = true;
for (int j = 0; j < f_nodes.size(); j++)
{
if (!equal(f_nodes[j].RealArray(g->mesh->CoordsTag())[2],minmax_z)) { done = false; break; }
}
if (done)
return faces[i];
}
return result_face;
}
void resolve_vertical_conflicts(grid* g)
{
Tag tag_owner = g->mesh->RedistributeTag();
for(Mesh::iteratorCell it = g->mesh->BeginCell(); it != g->mesh->EndCell(); it++)
{
Face uface = get_bottom_or_upper_face(g,it->getAsCell(),false);
if (uface.getAdjElements(CELL).size() <= 1)
{
Cell c = it->getAsCell();
int i = 0;
do
{
Face b_face = get_bottom_or_upper_face(g,c, true);
ElementArray<Element> adjs = b_face.getAdjElements(CELL);
if (adjs.size() <= 1) break;
if (adjs.size() != 2) TSNH
if (adjs[0]->getAsCell() != c) c = adjs[0]->getAsCell();
else c = adjs[1].getAsCell();
c.Integer(tag_owner) = it->Integer(tag_owner);
} while (i++ < 10000);
}
}
}
void old_correct_brothers(struct grid* g, int size, int rank, int type)
{
Tag tag_owner = g->mesh->RedistributeTag();
int* m = new int[size];
for(Mesh::iteratorCell it = g->mesh->BeginCell(); it != g->mesh->EndCell(); it++)
{
if (it->Integer(g->c_tags.level) == 0) continue;
for (int i = 0; i < size; i++) m[i] = 0;
int chld_num = it->Integer(g->c_tags.chld_num);
int center_node_num[8] = {6,5,1,2,7,4,0,3};
......@@ -1341,19 +1459,22 @@ void correct_brothers(struct grid* g, int size, int rank, int type)
ElementArray<Cell> children;
ElementArray<Element> adjs = center_node.getAdjElements(CELL);
//if (adjs.size() != 8) TSNH
if (adjs.size() != 8) { cout << rank << " " << adjs.size(); TSNH }
int my = 0;
int new_owner = 0;
for (ElementArray<Element>::iterator p = adjs.begin(); p != adjs.end(); p++)
{
if (p->getAsCell().Integer(tag_owner) == rank) my++;
int owner = p->getAsCell().Integer(tag_owner);
m[owner]++;
}
if (my < 4)
new_owner = (rank == 0 ? 1 : 0);
else
new_owner = (rank == 0 ? 0 : 1);
int max = 0;
for (int i = 0; i < size; i++)
if (max < m[i])
{
max = m[i];
new_owner = i;
}
for (ElementArray<Element>::iterator p = adjs.begin(); p != adjs.end(); p++)
{
......@@ -1361,8 +1482,56 @@ void correct_brothers(struct grid* g, int size, int rank, int type)
p->getAsCell().Integer(tag_owner) = new_owner;
}
}
resolve_vertical_conflicts(g);
}
void octree_refine(grid* g, int rank)
{
Tag tag_owner = g->mesh->RedistributeTag();
for(Mesh::iteratorCell it = g->mesh->BeginCell(); it != g->mesh->EndCell(); it++)
it->Integer(tag_owner) = rank;
}
typedef map<int,ElementArray<Cell> > cmap;
void correct_brothers(struct grid* g, int size, int rank, int type)
{
Tag tag_owner = g->mesh->RedistributeTag();
cmap cells_map;
for(Mesh::iteratorCell it = g->mesh->BeginCell(); it != g->mesh->EndCell(); it++)
{
(cells_map[it->Integer(g->c_tags.base_id)]).push_back(it->getAsCell());
}
int* m = new int[size];
for (cmap::iterator it = cells_map.begin(); it != cells_map.end(); it++)
{
for (int i = 0; i < size; i++) m[i] = 0;
for (ElementArray<Cell>::iterator p = it->second.begin(); p != it->second.end(); p++)
{
int owner = p->getAsCell().Integer(tag_owner);
m[owner]++;
}
int max = 0;
int new_owner = 0;
for (int i = 0; i < size; i++)
if (max < m[i])
{
max = m[i];
new_owner = i;
}
for (ElementArray<Cell>::iterator p = it->second.begin(); p != it->second.end(); p++)
p->getAsCell().Integer(tag_owner) = new_owner;
}
resolve_vertical_conflicts(g);
delete[] m;
}
void resolve_ghost(grid* g)
{
......@@ -1473,11 +1642,12 @@ void remove_border(grid* g)
// new_faces contains faces that should't be unite
// Unite 4 faces to one
/* OLD VERSION
ElementArray<Node> new_nodes;
ElementArray<Node> points_to_delete;
ElementArray<Element> l_edges = it->getAdjElements(EDGE);
for (ElementArray<Element>::iterator f = l_edges.begin(); f != l_edges.end(); f++)
{
if (f->getAsEdge().getBeg() == it->getAsNode()) points_to_delete.push_back(f->getAsEdge().getEnd());
......@@ -1506,7 +1676,23 @@ void remove_border(grid* g)
new_nodes.push_back(node->getAsNode());
}
}
*/
ElementArray<Node> new_nodes;
for (ElementArray<Element>::iterator f = adjs.begin(); f != adjs.end(); f++)
{
ElementArray<Node> nodes = f->getNodes();
for (ElementArray<Node>::iterator p = nodes.begin(); p != nodes.end(); p++)
if (p->getAsNode() != it->getAsNode())
{
int was = 0;
for (int i = 0; i < new_nodes.size(); i++)
if (new_nodes[i] == p->getAsNode()) was = 1;
if (was == 0)
new_nodes.push_back(p->getAsNode());
}
}
order_nodes_in_face_a(g,&new_nodes);
// Now we have 4 nodes for create new face
Face new_face = g->mesh->CreateFace(new_nodes).first;
......@@ -1660,15 +1846,17 @@ void resolve_connections(grid* g)
int yo = 0;
void gridAMR(struct grid * g, int action)
{
int rank = g->mesh->GetProcessorRank();
remove_border(g);
resolve_edges(g);
if (action == 0 || action == 1 || action == 3)
gridCoarse(g);
gridCoarse(g);
if (action == 0 || action == 2 || action == 3)
gridRefine(g);
gridRefine(g);
resolve_edges(g);
if (action == 3)
{
remove_border(g);
......
......@@ -44,11 +44,14 @@ struct Cell_Tags
/// Number of processor before redistribute
Tag proc;
Tag i;
Tag base_id;
};
/// Main object, contains all information about mesh
struct grid
{
int last_base_id;
/// Inmost mesh
Mesh* mesh;
Cell_Tags c_tags;
......@@ -61,6 +64,7 @@ struct grid
/// Often after redistribution brother cells splits to different processor
/// This function corrects this division. Puts all children to one processor
void old_correct_brothers(struct grid* g, int size, int rank, int type);
void correct_brothers(struct grid* g, int size, int rank, int type);
void pre_eval(struct grid* g, int size, int rank);
void init_mesh(struct grid* g);
......
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