diff --git a/Examples/Octree/main.cpp b/Examples/Octree/main.cpp index 84751b2ed29f75271709126356c9ece6b8a8dfac..da8246991e53a3932c32b96194530e576544b419 100644 --- a/Examples/Octree/main.cpp +++ b/Examples/Octree/main.cpp @@ -119,23 +119,10 @@ void rev_transformation(double xyz[3]) xyz[2] = (tmp[2] - 4010.0) / 10.0 + 0.5; } -/// Function provided to octgrid algorithm. -/// Defines that cell should be unite. Returns 1 for unite else returns 0. -int cell_should_unite(struct grid * g, Cell cell) -{ - const double r = base_radius; - int test = 1; - - double x = cell.RealArrayDF(g->c_tags.center)[0]; - double y = cell.RealArrayDF(g->c_tags.center)[1]; - - test &= (x-mx)*(x-mx)+(y-my)*(y-my) > r; - return test; -} /// Function provided to octgrid algorithm. /// Defines that cell should be split. Returns 1 for split else returns 0. -int cell_should_split(struct grid * g, Cell cell, int level) +int cell_should_split(struct grid * g, Cell cell) { double r = base_radius; @@ -151,7 +138,34 @@ int cell_should_split(struct grid * g, Cell cell, int level) for (int level = 2; level <= refine_depth; level++) { if ((x-mx)*(x-mx)+(y-my)*(y-my) < r*5*(level-1)) - if (c_level < refine_depth - level + 1) return 1; + if (c_level < refine_depth - level + 1) + { + return 1; + } + } + return 0; +} + + +/// Function provided to octgrid algorithm. +/// Defines that cell should be unite. Returns 1 for unite else returns 0. +int cell_should_unite(struct grid * g, Cell cell) +{ +// return !cell_should_split(g,cell); + const double r = base_radius; + + double x = cell.RealArrayDF(g->c_tags.center)[0]; + double y = cell.RealArrayDF(g->c_tags.center)[1]; + int c_level = cell.Integer(g->c_tags.level); + + if (c_level == refine_depth) + { + if ((x-mx)*(x-mx)+(y-my)*(y-my) > r) return 1; + } + else + { + double R = (refine_depth - c_level)*5*r; + if ((x-mx)*(x-mx)+(y-my)*(y-my) > R) return 1; } return 0; } diff --git a/Examples/Octree/main_light.cpp b/Examples/Octree/main_light.cpp index db715c4f941e31a768d0c8fedd943e64c418bf7e..6b08985260814aadd7425d17d178c56a3c03f992 100644 --- a/Examples/Octree/main_light.cpp +++ b/Examples/Octree/main_light.cpp @@ -96,7 +96,7 @@ int cell_should_unite(struct grid * g, Cell cell) /// Function provided to octgrid algorithm. /// Defines that cell should be split. Returns 1 for split else returns 0. -int cell_should_split(struct grid * g, Cell cell, int level) +int cell_should_split(struct grid * g, Cell cell) { double r = base_radius; diff --git a/Examples/Octree/main_test.cpp b/Examples/Octree/main_test.cpp index 0d280d5bcf735e83ad0a3efd65a0b4d95fc4ee3c..4d7e084fc7ef5428b6856e4fc13648b86d15205b 100644 --- a/Examples/Octree/main_test.cpp +++ b/Examples/Octree/main_test.cpp @@ -1,3 +1,5 @@ + +//make all && mpirun -np 1 ./Octree_test -log=3 -l=2 -n=50x50x2 #include "octgrid.h" #include #include "inmost.h" @@ -66,19 +68,28 @@ void pre_redistribute(int type) /// Defines that cell should be unite. Returns 1 for unite else returns 0. int cell_should_unite(struct grid * g, Cell cell) { +// return !cell_should_split(g,cell); const double r = base_radius; - int test = 1; double x = cell.RealArrayDF(g->c_tags.center)[0]; double y = cell.RealArrayDF(g->c_tags.center)[1]; + int c_level = cell.Integer(g->c_tags.level); - test &= (x-mx)*(x-mx)+(y-my)*(y-my) > r; - return test; + if (c_level == refine_depth) + { + if ((x-mx)*(x-mx)+(y-my)*(y-my) > r) return 1; + } + else + { + double R = (refine_depth - c_level)*5*r; + if ((x-mx)*(x-mx)+(y-my)*(y-my) > R) return 1; + } + return 0; } /// Function provided to octgrid algorithm. /// Defines that cell should be split. Returns 1 for split else returns 0. -int cell_should_split(struct grid * g, Cell cell, int level) +int cell_should_split(struct grid * g, Cell cell) { double r = base_radius; @@ -290,7 +301,7 @@ void parse_arguments(int argc, char** argv, int* n, double* R, int* L, int* log) } } -int iters_count = 5; +int iters_count = 6; int main(int argc, char ** argv) { @@ -325,6 +336,9 @@ int main(int argc, char ** argv) BARRIER double st = Timer(); double ct , tt; + double time_amr, time_red; + double a_amr = 0; + double a_red = 0; for (int iter = 0; iter < iters_count; iter++) { BARRIER @@ -333,22 +347,34 @@ int main(int argc, char ** argv) gridAMR(&thegrid,0); BARRIER tt = Timer(); - if (rank == 0) LOG(1, "AMR time = " << tt-ct); + time_amr = tt-ct; ct = tt; redistribute(&thegrid, 0); BARRIER tt = Timer(); - if (rank == 0) LOG(1, "Red time = " << tt-ct); + time_red = tt-ct; ct = tt; LOG(2, rank << ": iteration " << i << " complete. Cells: " << thegrid.mesh->NumberOfCells()) + if (iter > 0) + { + a_amr += time_amr; + a_red += time_red; + } + if (rank == 0) LOG(1, "AMR time = " << time_amr); + if (rank == 0) LOG(1, "Red time = " << time_red); + dump_to_vtk(&thegrid); + BARRIER + if (rank == 0) LOG(1, "==============="); i++; mx += h; } BARRIER - tt = Timer() - tt; + tt = Timer() - st; + if (rank == 0) cout << "time = " << tt << endl; + if (rank == 0) cout << "Average AMR time = " << a_amr/(iters_count - 1) << endl; + if (rank == 0) cout << "Average RED time = " << a_red/(iters_count - 1) << endl; if (rank == 0) cout << "time = " << tt << endl; - dump_to_vtk(&thegrid); // send_dump_command(); // send_quit_command(); diff --git a/Examples/Octree/octgrid.cpp b/Examples/Octree/octgrid.cpp index d0a1819bbc3b3c6aa17aca4287a6a9147571f78d..7c879acca3fdbc463e627f3f476c385a408b3aa8 100644 --- a/Examples/Octree/octgrid.cpp +++ b/Examples/Octree/octgrid.cpp @@ -3,6 +3,7 @@ #include #include #include +#include using namespace std; double epsilon = 0.0001; @@ -11,16 +12,110 @@ double epsilon = 0.0001; #define IS_GHOST_P(c) (c->GetStatus() == Element::Ghost) #define BARRIER MPI_Barrier(MPI_COMM_WORLD); +const int MAX_LEVEL = 10; void default_transformation(double xyz[3]) { (void) xyz; } int default_cell_should_unite(struct grid * g, int cell) { (void) g; (void) cell; return 0; } int default_cell_should_split(struct grid * g, int cell) { (void) g; (void) cell; return 0; } void mark_cell(struct grid* g, Cell cell); void cellCoarse(struct grid* g, Cell cell); void resolve_connections(grid* g); +void resolve_edges(grid* g, Cell* c = NULL); inline bool equal(double a, double b); void correct_brothers(struct grid* g, int size, int rank, int type); int global = 0; +double time_u; +double time_p; +double time_r1; +double time_r2; +double ttt; + + +class CR_cells +{ + ElementArray* d_cells[MAX_LEVEL]; + vector* d_valid[MAX_LEVEL]; + ElementArray d_cr_cells; // Array for refine and coarse + vector d_cr_valids; + bool d_empty; + grid* g; + + public: + bool active; + CR_cells() + { + d_empty = true; + active = false; + for (int i = 0; i < MAX_LEVEL; i++) d_cells[i] = NULL; + for (int i = 0; i < MAX_LEVEL; i++) d_valid[i] = NULL; + } + void push(Cell c) + { + int l = c.Integer(g->c_tags.level); + if (d_cells[l] == NULL) d_cells[l] = new ElementArray; + if (d_valid[l] == NULL) d_valid[l] = new vector; + + d_cells[l]->push_back(c); + d_valid[l]->push_back(true); + c.Integer(g->c_tags.i) = d_cells[l]->size() - 1; + +/* + d_cr_cells.push_back(c); + d_cr_valids.push_back(true); + */ + d_empty = false; + } + void set_grid(grid* _g) { g = _g; } + +// bool check(int i) { return i >= d_cr_valids.size() ? false : d_cr_valids[i]; } + bool check(int l, int i) { return d_valid[l] == NULL ? false : ( i >= d_valid[l]->size() ? false : (*d_valid[l])[i]); } + + // void dis(int i) { d_cr_valids[i] = false; } + void dis(int l, int i) { (*d_valid[l])[i] = false; } + + // int size() { return d_cr_cells.size(); } + int size(int l) { return d_cells[l] == NULL ? 0 : d_cells[l]->size(); } + + // Cell operator[](int i) { return d_cr_cells[i]; } + Cell get(int l, int i) + { + if (d_cells[l] == NULL) TSNH; + return (*(d_cells[l]))[i]; + } + + void clear() + { + d_cr_cells.clear(); d_cr_valids.clear(); + for (int i = 0; i < MAX_LEVEL; i++) + { + if (d_cells[i] != NULL) delete d_cells[i]; + if (d_valid[i] != NULL) delete d_valid[i]; + d_cells[i] = NULL; + d_valid[i] = NULL; + } + d_empty = true; + } + + // void set(Cell c, int i) { d_cr_cells[i] = c; } + void set1(Cell c, int i) + { + int l = c.Integer(g->c_tags.level); + (*d_cells[l])[i] = c; + c.Integer(g->c_tags.i) = i; + } + + bool empty() { return d_empty; } + // bool check_i(int i ) { return (i >= 0 && i < d_cr_cells.size()); } + bool check_i(int l, int i ) { return d_cells[l] == NULL ? false : (i >= 0 && i < d_cells[l]->size()); } + + void print_size() + { + for (int i = 0; i < MAX_LEVEL; i++) + if (d_cells[i] != NULL) cout << i << ": " << d_cells[i]->size() << endl; + } +}; +CR_cells cr_cells; + /// Dump mesh to vtk file in folder "grids" void dump_to_vtk(grid* g) { @@ -335,6 +430,7 @@ void gridInit(struct grid * g, int n[3]) int i,j,k,l,m, c[3], v, r, env[1<mesh = new Mesh(); g->mesh->SetDimensions(3); @@ -558,11 +654,27 @@ Cell recreate_cell(struct grid* g, Cell cell, ElementArray faces) res.Integer(g->c_tags.floor) = cell.Integer(g->c_tags.floor); res.Integer(g->c_tags.chld_num) = cell.Integer(g->c_tags.chld_num); + int l = res.Integer(g->c_tags.level) = cell.Integer(g->c_tags.level); + int cr_i = res.Integer(g->c_tags.i) = cell.Integer(g->c_tags.i); + +/* + if (cr_cells.check_i(l,cr_i)) + if (cr_cells.get(l.cr_i) == cell) + cr_cells.set1(res,cr_i); +*/ + if (cr_cells.active && l > 0) + if (cr_cells.check_i(l,cr_i)) + { + if (cr_cells.get(l,cr_i) == cell) + cr_cells.set1(res,cr_i); + else + TSNH + } + for (int i = 0; i < 8; i++) res.IntegerArrayDF(g->c_tags.par_chld_nums)[i] = cell.IntegerArrayDF(g->c_tags.par_chld_nums)[i]; - 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); @@ -1040,8 +1152,15 @@ Cell cell_unite(struct grid* g, ElementArray children, Node center) order_nodes_in_face(g, &def_nodes); + // Destroy cells for (ElementArray::iterator p = children.begin(); p != children.end(); p++) + { + int cr_i = p->Integer(g->c_tags.i); + int l = p->Integer(g->c_tags.level); + //cr_cells.dis(cr_i); + cr_cells.dis(l,cr_i); g->mesh->Destroy(p->getAsCell()); + } g->mesh->Destroy(center); int x_nodes[6] = { 4,15,21,10,12,13 }; @@ -1062,7 +1181,7 @@ Cell cell_unite(struct grid* g, ElementArray children, Node center) if ((new_faces.size() - m) != 4 ) TSNH } - // This case means that adjacent cell in on cell. + // This case means that adjacent cell in one cell. if (adjs_count == 1) { ElementArray face_nodes; @@ -1121,7 +1240,6 @@ Cell cell_unite(struct grid* g, ElementArray children, Node center) } Cell c = g->mesh->CreateCell(new_faces).first; - g->rev_transformation(xyz); c.RealArrayDF(g->c_tags.center)[0] = xyz[0]; c.RealArrayDF(g->c_tags.center)[1] = xyz[1]; @@ -1153,6 +1271,20 @@ Cell cell_unite(struct grid* g, ElementArray children, Node center) break; } } + +/* + int cr_i = cr_cells.size(); + c.Integer(g->c_tags.i) = cr_i; + cr_cells.push(c); +*/ + + if (level - 1 > 0) + { + int cr_i = cr_cells.size(level - 1); + c.Integer(g->c_tags.i) = cr_i; + cr_cells.push(c); + } + return c; } @@ -1161,7 +1293,7 @@ Cell cell_unite(struct grid* g, ElementArray children, Node center) bool cell_check_mark(struct grid* g, Cell cell) { if (cell->Integer(g->c_tags.to_split) == 1) return false; - if( g->cell_should_split(g,cell,2) ) + if( g->cell_should_split(g,cell) ) { mark_cell(g,cell); return true; @@ -1284,7 +1416,6 @@ void cellCoarse(struct grid* g, Cell cell) if (adjs.size() != 8) { cout << adjs.size(); TSNH } for (ElementArray::iterator adj = adjs.begin(); adj != adjs.end(); adj++) { - ElementArray local_nodes = adj->getAsCell().getNodes(); for (ElementArray::iterator p = local_nodes.begin(); p != local_nodes.end(); p++) { @@ -1292,7 +1423,19 @@ void cellCoarse(struct grid* g, Cell cell) for (ElementArray::iterator q = local_adjs.begin(); q != local_adjs.end(); q++) if (q->getAsCell().Integer(g->c_tags.level) > adj->getAsCell().Integer(g->c_tags.level)) { + int cr_i = cell.Integer(g->c_tags.i); + int l = cell.Integer(g->c_tags.level); cellCoarse(g,q->getAsCell()); + //if (cr_cells.check(i)) + //{ + // cr_cells.push(cr_cells[i]); + // cr_cells.dis(i); + // } + if (cr_cells.check(l,i)) + { + cr_cells.push(cr_cells.get(l,i)); + cr_cells.dis(l,i); + } return; } } @@ -1300,7 +1443,14 @@ void cellCoarse(struct grid* g, Cell cell) children.push_back(adj->getAsCell()); } + time_p += Timer() - ttt; + ttt = Timer(); Cell c = cell_unite(g,children, center_node); + time_u += Timer() - ttt; + ttt = Timer(); + resolve_edges(g,&c); + time_r1 += Timer() - ttt; + ttt = Timer(); return; } @@ -1391,21 +1541,31 @@ void unite_2_edges(struct grid* g, Node node) g->mesh->Destroy(node); } -void resolve_edges(grid* g) +void resolve_edges(grid* g, Cell* c) { - for(Mesh::iteratorEdge it = g->mesh->BeginEdge(); it != g->mesh->EndEdge(); it++) - { - ElementArray adjs = it->getAdjElements(FACE); - } - - for(Mesh::iteratorNode it = g->mesh->BeginNode(); it != g->mesh->EndNode(); it++) - { - ElementArray edges = it->getEdges(); - if (edges.size() == 2) + if (c) + { + ElementArray nodes = c->getNodes(); + for (int i = 0; i < nodes.size(); i++) + if (nodes[i].getEdges().size() == 2) { - unite_2_edges(g,it->getAsNode()); + time_r1 += Timer() - ttt; + ttt = Timer(); + unite_2_edges(g,nodes[i]); + time_r2 += Timer() - ttt; + ttt = Timer(); } + return; + } + + for(Mesh::iteratorNode it = g->mesh->BeginNode(); it != g->mesh->EndNode(); it++) + { + ElementArray edges = it->getEdges(); + if (edges.size() == 2) + { + unite_2_edges(g,it->getAsNode()); } + } } void gridCoarse(struct grid * g) @@ -1413,6 +1573,66 @@ void gridCoarse(struct grid * g) int rank = g->mesh->GetProcessorRank(); int m = 0; bool ready; + cr_cells.clear(); + bool empty = true; + cr_cells.active = true; + double tt = Timer(); + int i = 0; + for(Mesh::iteratorCell it = g->mesh->BeginCell(); it != g->mesh->EndCell(); it++) + { + if (it->Integer(g->c_tags.level) == 0) continue; + // if (g->cell_should_unite(g,it->getAsCell())) + { + empty = false; + it->Integer(g->c_tags.i) = i; + cr_cells.push(it->getAsCell()); + i++; + } + // else + // it->Integer(g->c_tags.i) = -1; + } + cr_cells.print_size(); + cout << "N time: " << Timer() - tt << endl; + if (empty) + { + cr_cells.push(g->mesh->BeginCell()->getAsCell()); + return; + } + + time_u = 0; + time_p = 0; + time_r1 = 0; + time_r2 = 0; + ttt = Timer(); + for (int l = MAX_LEVEL-1; l >= 1; l--) + { + for (int i = 0; i < cr_cells.size(l); i++) + { + if (cr_cells.check(l,i) == false) continue; + if (!cr_cells.get(l,i).isValid()) continue; + if (cr_cells.get(l,i).GetStatus() == Element::Ghost) continue; + // cell_check_coarse(g,cr_cells.get(l,i)); + Cell c = cr_cells.get(l,i); + if (g->cell_should_unite(g,c)) + cellCoarse(g,c); + + } + } + /* + for (int i = 0; i < cr_cells.size(); i++) + { + if (cr_cells.check(i) == false) continue; + if (!cr_cells[i].isValid()) continue; + if (cr_cells[i].GetStatus() == Element::Ghost) continue; + cell_check_coarse(g,cr_cells[i]); + } + */ + time_p += Timer() - ttt; + cout << "Time Prep: " << time_p << endl; + cout << "Time Unit: " << time_u << endl; + cout << "Time Res1: " << time_r1 << endl; + cout << "Time Res2: " << time_r2 << endl; + /* while (1) { ready = true; @@ -1426,9 +1646,12 @@ void gridCoarse(struct grid * g) } } - resolve_edges(g); + //resolve_edges(g); if (ready) break; } + */ + + cr_cells.active = false; } void print_redist_tag(struct grid* g, int rank) @@ -1855,7 +2078,7 @@ bool split_if_needed(grid* g, Cell c, Face f) c.RealArrayDF(g->c_tags.center)[1] = fc[1]; c.RealArrayDF(g->c_tags.center)[2] = fc[2]; - int res = g->cell_should_split(g, c , 2); + int res = g->cell_should_split(g, c ); c.RealArrayDF(g->c_tags.center)[0] = cx; c.RealArrayDF(g->c_tags.center)[1] = cy; c.RealArrayDF(g->c_tags.center)[2] = cz; @@ -1915,48 +2138,74 @@ void resolve_connections(grid* g) } } -double tt; -void start_time() -{ - BARRIER - double tt = Timer(); -} - -void print_time(string what, int rank) +class IterTime { - BARRIER - double tt1 = Timer(); - if (rank == 0) cout << setw(15) << what << " time = " << tt1 - tt << endl; - tt = tt1; -} + vector time; + vector what; + double tt; + double all; + + public: + void start() + { + BARRIER + tt = Timer(); + } + void add(string _what) + { + double tt1 = Timer(); + double ct = tt1 - tt; + all += ct; + time.push_back(ct); + what.push_back(_what); + tt = tt1; + } + void print_time(int rank) + { + if (rank > 0) return; + for (int i = 0; i < time.size(); i++) + { + double per = (time[i]/all)*100; + cout << setw(15) << what[i] << " time = " << setw(10) << time[i] << " : "; + if (per > 70) + cout << "\033[1;31m" << setprecision(2) << setw(5) << (time[i]/all)*100 << " %" << "\033[0m" << endl; + else + cout << setprecision(2) << setw(5) << (time[i]/all)*100 << " %" << endl; + } + } +}; int yo = 0; void gridAMR(struct grid * g, int action) { int rank = g->mesh->GetProcessorRank(); // Get the rank of the current process - start_time(); + IterTime time; + time.start(); + remove_border(g); - print_time("Remove border",rank); + time.add("Remove border"); resolve_edges(g); - print_time("Resolve edges",rank); + time.add("Resolve edges"); if (action == 0 || action == 1 || action == 3) gridCoarse(g); - print_time("Coarse",rank); + time.add("Coarse"); - if (action == 0 || action == 2 || action == 3) + //if (action == 0 || action == 2 || action == 3) + if (action != 3) gridRefine(g); - print_time("Refine",rank); + time.add("Refine"); resolve_edges(g); - print_time("Resolve edges",rank); + time.add("Resolve edges"); - if (action == 3) + if (action == 4) { remove_border(g); resolve_edges(g); resolve_connections(g); } g->mesh->ResolveShared(); // Resolve duplicate nodes + time.print_time(rank); } diff --git a/Examples/Octree/octgrid.h b/Examples/Octree/octgrid.h index c730c8988279342da975e54093ae9a7b3c725593..487f3d068f131225f04062f1c4c8ba717a876e7b 100644 --- a/Examples/Octree/octgrid.h +++ b/Examples/Octree/octgrid.h @@ -58,7 +58,7 @@ struct grid void (*transformation)(double xyz[3]); void (*rev_transformation)(double xyz[3]); - int (*cell_should_split)(struct grid * g, Cell cell, int level); + int (*cell_should_split)(struct grid * g, Cell cell); int (*cell_should_unite)(struct grid * g, Cell cell); };