Added test

Examples/Octree/CMakeLists.txt

@@ -10,9 +10,15 @@ set(SOURCE_LIGHT main_light.cpp
 	   octgrid.cpp
 	   octgrid.h)
 
+set(SOURCE_TEST main_test.cpp
+	   octgrid.cpp
+	   octgrid.h)
+
+
 add_executable(Octree ${SOURCE})
 add_executable(Octree_light ${SOURCE_LIGHT})
-set(TARGETS Octree Octree_light)
+add_executable(Octree_test ${SOURCE_TEST})
+set(TARGETS Octree Octree_light Octree_test)
 
 foreach(vtarget ${TARGETS})
 target_link_libraries(${vtarget} inmost)

Examples/Octree/main.cpp

@@ -26,7 +26,7 @@ int cur_cell = 0; bool draw_all_cells = true;
 int cur_face = 0; bool draw_all_faces = true;
 int cur_edge = 0; bool draw_all_edges = true;
 
-int from_file = 0;
+int from_file = 1;
 int redist_after_amr = 0;
 
 // Variables for refine and coarse
@@ -75,22 +75,7 @@ double** sub_edges_nodes;
 int*     sub_edges_nodes_n; 
 char**    sub_edges_ghost;
 
-/// Dump mesh to vtk file in folder "grids"
-void dump_to_vtk()
-{
-	//thegrid.mesh->ResolveShared(); // Resolve duplicate nodes
-	//thegrid.mesh->ExchangeGhost(2,NODE); // Construct Ghost cells in 2 layers connected via nodes
-
-    std::stringstream filename;
-    filename << "grids/grid_";
-    filename << size;
-    if( size == 1 )
-        filename << ".vtk";
-    else
-        filename << ".pvtk";
-    thegrid.mesh->Save(filename.str());
-	cout << "Process " << ::rank << ": dumped mesh to file" << endl;
-}
+
 
 void print_all_cells(grid* g)
 {
@@ -526,54 +511,7 @@ void prepare_to_correct_brothers()
     thegrid.mesh->AssignGlobalID(CELL | EDGE | FACE | NODE);
 }
 
-/// 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 = 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);
-    
-    try
-    {
-        part->Evaluate();
-    }
-    catch (INMOST::ErrorType er)
-    {
-        cout << "Exception: " << er << endl;
-    }
-    catch(...)
-    {
-    }
-    delete part;
 
-    thegrid.mesh->RemoveGhost();
-	correct_brothers(&thegrid,size,::rank, 2);
-
-    try
-    {
-        thegrid.mesh->Redistribute(); 
-    }
-    catch (INMOST::ErrorType er)
-    {
-        cout << "Exception: " << er << endl;
-    }
-    catch(...)
-    {
-    }
-
-
-    thegrid.mesh->RemoveGhost();
-    thegrid.mesh->ReorderEmpty(CELL|FACE|EDGE|NODE);
-    //thegrid.mesh->AssignGlobalID(CELL | EDGE | FACE | NODE);
-	LOG(2,"Process " << ::rank << ": redistribute completed")
-}
 
 /// Prepare to redistribute. Master sends special command to slaves which means redistribute
 void redistribute_command()
@@ -596,7 +534,7 @@ void redistribute_command()
         MPI_Isend(buff[i], 2, MPI_CHAR, i, 0, INMOST_MPI_COMM_WORLD, req + i);
     }       
     
-    redistribute(type);
+    redistribute(&thegrid, type);
     counter++;
 
 }
@@ -727,7 +665,7 @@ void keyboard(unsigned char key, int x, int y)
     if( key == 'f' ) 
 	{
 		if (::rank == 0) { // Send dump command to other process
-            dump_to_vtk();
+            dump_to_vtk(&thegrid);
             char buff[10][2];
             MPI_Request req[10];
             for (int i = 1; i < size; i++)
@@ -812,7 +750,7 @@ void NotMainProcess()
         {
             char type_C = buff[1];
             int  type = type_C - '0';
-            redistribute(type);
+            redistribute(&thegrid,type);
         }
         if (buff[0] == 'u') // Need remove ghosts
         {
@@ -828,7 +766,7 @@ void NotMainProcess()
         }
         if (buff[0] == 'f') // Dump mesh to file
         {
-            dump_to_vtk();
+            dump_to_vtk(&thegrid);
         }
         if (buff[0] == 'r') // Master wants to redraw mesh. Now this procces need send own grid to master.
         {
@@ -887,7 +825,7 @@ void NotMainProcess()
     }
 }
 
-void parse_arguments(int argc, char** argv, int* n, double* R, int* L)
+void parse_arguments(int argc, char** argv, int* n, double* R, int* L, string& file_name)
 {
   if (argc < 2) return;
 
@@ -930,6 +868,10 @@ void parse_arguments(int argc, char** argv, int* n, double* R, int* L)
     {
 	    *L = atoi(str2.c_str());
     }
+    else if (str1 == "-f")
+    {
+        file_name = str2;
+    }
     else
     {
       cout << "Invalid command: " << str1 << endl;
@@ -944,6 +886,7 @@ void print_help()
   cout << "   -n=10x10x1     - grid size" << endl;
   cout << "   -r=0.01        - refine radius" << endl;
   cout << "   -l=2           - refine level"  << endl;
+  cout << "   -f=grid_3.pvtk - takes grid from file" << endl;
   cout << endl;
   cout << "Hotkeys:" << endl;
   cout << "   Space - refine grid around mouse cursor" << endl;
@@ -967,10 +910,21 @@ int main(int argc, char ** argv)
 	MPI_Comm_rank(MPI_COMM_WORLD, &::rank);
 
     if (::rank == 0) print_help();
-    parse_arguments(argc, argv, n, &base_radius,&refine_depth);
+    string file_name;
+    parse_arguments(argc, argv, n, &base_radius,&refine_depth, file_name);
     all_cells_count = n[0]*n[1]*n[2] * 2;
     
+    if (file_name.empty()) {
         gridInit(&thegrid,n);
+    } else {   
+        thegrid.mesh = new Mesh(); // Create an empty mesh
+        thegrid.mesh->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
+        file_name = "grids/" + file_name;
+        thegrid.mesh->Load(file_name.c_str()); // Load input mesh
+        //	thegrid.mesh->Load("dump.pvtk"); // Load input mesh
+        init_mesh(&thegrid); 
+    }
+
     Partitioner::Initialize(&argc,&argv);
 
 	::size = thegrid.mesh->GetProcessorsNumber();

Examples/Octree/main_test.cpp

+#include "octgrid.h"
+#include <math.h>
+#include "inmost.h"
+#include <iomanip>
+#include <iostream>
+
+#define LOG(level,msg)  { if (log_level >= level) cout << msg << endl; }
+#define BARRIER MPI_Barrier(MPI_COMM_WORLD);
+
+using namespace std;
+
+const int MAX_PROCESSORS_COUNT = 10;
+
+struct grid thegrid;
+int refine_depth = 2;
+double  mx = 0,  my = 0;
+double rmx = 0, rmy = 0;
+double base_radius = 0.01;
+int action = 0;
+int log_level = 0;
+int all_cells_count;
+int rank;
+int size;
+int counter = 0; // Count number of steps
+
+/// Function provided to octgrid algorithm. Defines transformation from grid to grid for draw.
+void transformation(double xyz[3]) 
+{
+    double tmp[3];
+	tmp[0] = xyz[0];
+	tmp[1] = xyz[1];
+	tmp[2] = xyz[2];
+	tmp[0] -= 0.5;
+	tmp[0] *= 100.0;
+	tmp[1] -= 0.5;
+	tmp[1] *= 100.0;
+	xyz[0] =  tmp[0];
+	xyz[1] =  tmp[1];
+	//xyz[2] = 4010.0 + 10.0 * (tmp[2]-0.5);
+	xyz[2] = 4010.0 + 10.0 * (tmp[2]-0.5);
+}
+
+/// Function provided to octgrid algorithm. Defines transformation from grid for draw to grid.
+void rev_transformation(double xyz[3]) 
+{
+    double tmp[3];
+	tmp[0] = xyz[0];
+	tmp[1] = xyz[1];
+	tmp[2] = xyz[2];
+	tmp[0] /= 100.0;
+	tmp[0] += 0.5;
+	tmp[1] /= 100.0;
+	tmp[1] += 0.5;
+	xyz[0] =  tmp[0];
+	xyz[1] =  tmp[1];
+	xyz[2] = 4010.0 + 10.0 * (tmp[2]-0.5);
+	xyz[2] = (tmp[2] - 4010.0) / 10.0 + 0.5;
+}
+
+void pre_redistribute(int type)
+{
+    redistribute(&thegrid, type);
+}
+
+/// 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)
+{
+	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 ((x-mx)*(x-mx)+(y-my)*(y-my) < r)
+    {
+        if (c_level < refine_depth) return 1;
+    }
+
+    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;
+    }
+    return 0;
+}
+
+/// Main loop for not main process
+void NotMainProcess()
+{
+
+    char buff[256];
+    MPI_Status status;
+    int length;
+    int iteration = 0;
+
+    while (1) 
+    {
+        MPI_Recv(buff, 256, MPI_CHAR, 0,0, INMOST_MPI_COMM_WORLD, &status);
+
+		LOG(2, "Process " << ::rank << ": received message '" << buff[0] << "'")
+        
+        if (buff[0] == 'q') // Need to refine, mouse coordinates come
+		{
+			break;
+
+		}
+        if (buff[0] == 'm') // Need to refine, mouse coordinates come
+        {
+            mx = *((double*)(buff + 1));
+            my = *((double*)(buff + 1 + sizeof(double)));
+            int action = *((int*)(buff + 1 + 2*sizeof(double)));
+
+            gridAMR(&thegrid, action); 
+            //thegrid.mesh->AssignGlobalID(CELL | EDGE | FACE | NODE);
+        }
+        if (buff[0] == 'x') // Need to redistribute 
+        {
+            char type_C = buff[1];
+            int  type = type_C - '0';
+            //redistribute(&thegrid,type);
+            pre_redistribute(type);
+            cout << rank << ": iteration " << iteration++ << " complete. Cells: " << thegrid.mesh->NumberOfCells() << endl;
+        }
+        if (buff[0] == 'u') // Need remove ghosts
+        {
+            thegrid.mesh->RemoveGhost();
+        }
+        if (buff[0] == 'c') // Need command
+        {
+            command(&thegrid);
+        }
+        if (buff[0] == 'f') // Dump mesh to file
+        {
+            dump_to_vtk(&thegrid);
+        }
+   }
+}
+
+void send_dump_command()
+{
+	char buff[10][2];
+	MPI_Request req[10];
+	for (int i = 1; i < size; i++)
+	{
+		buff[i][0] = 'f'; // Special key, means dump file
+		MPI_Isend(buff[i], 1, MPI_CHAR, i, 0, INMOST_MPI_COMM_WORLD, req + i);
+	}
+}
+
+void send_quit_command()
+{
+	char buff[10][2];
+	MPI_Request req[10];
+	for (int i = 1; i < size; i++)
+	{
+		buff[i][0] = 'q'; // Special key, means dump file
+		MPI_Isend(buff[i], 1, MPI_CHAR, i, 0, INMOST_MPI_COMM_WORLD, req + i);
+	}
+}
+
+/// Send mpuse coordinates to slaves. Slaves receive the message 'm', execute "refine" with received coordinates.
+void send_coordinates_to_slaves(int action)
+{
+    char buff[MAX_PROCESSORS_COUNT][10];
+    MPI_Request req[MAX_PROCESSORS_COUNT];
+
+    for (int i = 1; i < size; i++)
+    {
+        LOG(2,"Main process: send coordinates to " << i)
+        buff[i][0] = 'm'; // Special key, means refine 
+        *((double*)(buff[i] + 1)) = mx;
+        *((double*)(buff[i] + 1 + sizeof(double))) = my;
+        *((int*)(buff[i] + 1 + 2*sizeof(double))) = action;
+        MPI_Isend(buff[i], 1 + 2 * sizeof(double) + sizeof(int), MPI_CHAR, i, 0, INMOST_MPI_COMM_WORLD, req + i);
+    }
+}
+
+
+
+/// Prepare to redistribute. Master sends special command to slaves which means redistribute
+void redistribute_command()
+{
+    char buff[MAX_PROCESSORS_COUNT][10];
+    MPI_Request req[MAX_PROCESSORS_COUNT];
+    int type = 2;
+    
+    if (counter == 5)
+    {
+        counter = 0;
+        type = 1;
+    }
+   
+    for (int i = 1; i < size; i++)
+    {
+		LOG(3,"Master: send redistribute command to slave " << ::rank)
+        buff[i][0] = 'x'; // Special key, means redistribute
+        buff[i][1] = type + '0'; // Special key, means redistribute
+        MPI_Isend(buff[i], 2, MPI_CHAR, i, 0, INMOST_MPI_COMM_WORLD, req + i);
+    }       
+    
+//    redistribute(&thegrid, type);
+    pre_redistribute(type);
+    counter++;
+
+}
+
+
+void print_help()
+{
+  cout << "Example of Octree refine on redistributed grid" << endl;
+  cout << "Command arguments:" << endl;
+  cout << "   -n=10x10x1 - grid size" << endl;
+  cout << "   -r=0.01    - refine radius" << endl;
+  cout << "   -l=2       - refine level"  << endl;
+  cout << "   -log=1     - log level"  << 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;
+}
+
+void parse_arguments(int argc, char** argv, int* n, double* R, int* L, int* log)
+{
+  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 if (str1 == "-l")
+    {
+	    *L = atoi(str2.c_str());
+    }
+    else if (str1 == "-log")
+    {
+	    *log = atoi(str2.c_str());
+    }
+    else
+    {
+      cout << "Invalid command: " << str1 << endl;
+    }
+  }
+}
+
+int iters_count = 5;
+
+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;
+	thegrid.cell_should_unite = cell_should_unite;
+    Mesh::Initialize(&argc,&argv);
+	MPI_Comm_rank(MPI_COMM_WORLD, &::rank);
+
+    if (::rank == 0) print_help();
+    parse_arguments(argc, argv, n, &base_radius,&refine_depth,&log_level);
+    all_cells_count = n[0]*n[1]*n[2] * 2;
+
+    gridInit(&thegrid,n);
+    Partitioner::Initialize(&argc,&argv);
+
+	::size = thegrid.mesh->GetProcessorsNumber();
+	::rank = thegrid.mesh->GetProcessorRank();
+
+    //dump_to_vtk();
+	if (rank == 0) cout << "Test start" << endl;
+
+    {
+		mx = 0.1;
+		my = 0.5;
+        double h = 0.8 / iters_count;
+        int i = 0;
+   		BARRIER
+		double st = Timer();
+        double ct , tt;
+        for (int iter = 0; iter < iters_count; iter++)
+        {
+    		BARRIER
+            ct = Timer();
+            if (rank == 0) LOG(1, "Iteration: " << i)
+            gridAMR(&thegrid,0);
+    		BARRIER
+            tt = Timer();
+    		if (rank == 0) LOG(1, "AMR time = " << tt-ct);
+            ct = tt;
+            redistribute(&thegrid, 0);
+    		BARRIER
+            tt = Timer();
+    		if (rank == 0) LOG(1, "Red time = " << tt-ct);
+            ct = tt;
+            LOG(2, rank << ": iteration " << i << " complete. Cells: " << thegrid.mesh->NumberOfCells())
+			i++;
+            mx += h;
+        }
+    	BARRIER
+		tt = Timer() - tt;
+		if (rank == 0) cout << "time = " << tt << endl;
+
+		dump_to_vtk(&thegrid);
+//		send_dump_command();
+//		send_quit_command();
+
+    }
+
+	Mesh::Finalize();
+}

Examples/Octree/octgrid.cpp

@@ -7,8 +7,8 @@ using namespace std;
 
 double epsilon = 0.0001;
 
-#define IS_GHOST(c) c.GetStatus() == Element::Ghost
-#define IS_GHOST_P(c) c->GetStatus() == Element::Ghost
+#define IS_GHOST(c) (c.GetStatus() == Element::Ghost)
+#define IS_GHOST_P(c) (c->GetStatus() == Element::Ghost)
 #define BARRIER MPI_Barrier(MPI_COMM_WORLD);
 
 void default_transformation(double xyz[3]) { (void) xyz; }
@@ -18,17 +18,78 @@ 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);
+void correct_brothers(struct grid* g, int size, int rank, int type);
 int global = 0;
 
+/// Dump mesh to vtk file in folder "grids"
 void dump_to_vtk(grid* g)
 {
 	//thegrid.mesh->ResolveShared(); // Resolve duplicate nodes
 	//thegrid.mesh->ExchangeGhost(2,NODE); // Construct Ghost cells in 2 layers connected via nodes
+    int rank = g->mesh->GetProcessorRank(); // Get the rank of the current process
+	int size = g->mesh->GetProcessorsNumber();
 
     std::stringstream filename;
-    filename << "grids/grid";
+    filename << "grids/grid_";
+    filename << size;
+    if( size == 1 )
+        filename << ".vtk";
+    else
         filename << ".pvtk";
     g->mesh->Save(filename.str());
+	cout << "Process " << rank << ": dumped mesh to file" << endl;
+}
+
+
+/// Redistribute grid by  partitioner
+void redistribute(grid* g, int type)
+{
+    g->mesh->RemoveGhost();
+    int rank = g->mesh->GetProcessorRank(); // Get the rank of the current process
+	int size = g->mesh->GetProcessorsNumber();
+
+	//LOG(2,"Process " << rank << ": redistribute. Cells: " << g->mesh->NumberOfCells())
+    Partitioner * part = new Partitioner(g->mesh);
+    
+    // Specify the partitioner
+    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);
+    
+    try
+    {
+        part->Evaluate();
+    }
+    catch (INMOST::ErrorType er)
+    {
+        cout << "Exception: " << er << endl;
+    }
+    catch(...)
+    {
+    }
+    delete