main.cpp 9.26 KB
Newer Older
1
2
3
4
5
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

Igor Konshin's avatar
Igor Konshin committed
6
#include "../../inmost.h"
7
using namespace INMOST;
Igor Konshin's avatar
Igor Konshin committed
8

9
10
11
12
13
14
15
16
17
#define MESH_SIZE 64

#define V_ID(x, y, z) ((x-localstart[0])*(localsize[1]+1)*(localsize[2]+1) + (y-localstart[1])*(localsize[2]+1) + (z-localstart[2]))

Mesh * ParallelCubeGenerator(INMOST_MPI_Comm comm, int nx, int ny, int nz)
{
	int procs_per_axis[3] = {1,1,1};
	int sizes[3] = {nx,ny,nz};
	int rank,size;
Igor Konshin's avatar
Igor Konshin committed
18
19
20
	Mesh * m = new Mesh(); // Create a mesh to be constructed

	m->SetCommunicator(comm); // Set the MPI communicator, usually MPI_COMM_WORLD
21
22

#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
23
	//MPI::COMM_WORLD.Set_errhandler(MPI::ERRORS_THROW_EXCEPTIONS);
24
#endif
Igor Konshin's avatar
Igor Konshin committed
25
26
27
28
29

	rank = m->GetProcessorRank(); // Get the rank of the current process
	size = m->GetProcessorsNumber(); // Get the number of processors used in communicator comm

	// Compute the configuration of processes connection
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
	{
		int divsize = size;
		std::vector<int> divs;
		while( divsize > 1 )
		{
			for(int k = 2; k <= divsize; k++)
				if( divsize % k == 0 )
				{
					divs.push_back(k);
					divsize /= k;
					break;
				}
		}
		int elements_per_procs[3] = {sizes[0],sizes[1],sizes[2]};
		for(std::vector<int>::reverse_iterator it = divs.rbegin(); it != divs.rend(); it++)
		{
			int * max = std::max_element(elements_per_procs+0,elements_per_procs+3);
			procs_per_axis[max-elements_per_procs] *= *it;
			(*max) /= *it;
		}
	}
Igor Konshin's avatar
Igor Konshin committed
51

52
	//rank = proc_coords[2] * procs_per_axis[0] *procs_per_axis[1] + proc_coords[1] * procs_per_axis[0] + proc_coords[0];
Igor Konshin's avatar
Igor Konshin committed
53
54
	int proc_coords[3] = {rank % procs_per_axis[0] , rank / procs_per_axis[0] % procs_per_axis[1], rank / (procs_per_axis[0] *procs_per_axis[1]) };

55
56
	//???
	int localsize[3], localstart[3], localend[3];
Igor Konshin's avatar
Igor Konshin committed
57
	int avgsize[3] =
58
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
59
60
61
		(int)ceil((double)sizes[0]/procs_per_axis[0]),
		(int)ceil((double)sizes[1]/procs_per_axis[1]),
		(int)ceil((double)sizes[2]/procs_per_axis[2])
62
	};
Igor Konshin's avatar
Igor Konshin committed
63

64
65
66
67
68
69
70
71
	for(int j = 0; j < 3; j++)
	{
		localstart[j] = avgsize[j] * proc_coords[j];
		if( proc_coords[j] == procs_per_axis[j] - 1 )
			localsize[j] = sizes[j] - avgsize[j] * (procs_per_axis[j]-1);
		else localsize[j] = avgsize[j];
		localend[j] = localstart[j] + localsize[j];
	}
Igor Konshin's avatar
Igor Konshin committed
72

73
74
	std::vector<Node *> newverts;
	newverts.reserve(localsize[0]*localsize[1]*localsize[2]);
Igor Konshin's avatar
Igor Konshin committed
75

76
77
78
79
80
81
82
83
	for(int i = localstart[0]; i <= localend[0]; i++)
		for(int j = localstart[1]; j <= localend[1]; j++)
			for(int k = localstart[2]; k <= localend[2]; k++)
			{
				Storage::real xyz[3];
				xyz[0] = i * 1.0 / (sizes[0]);
				xyz[1] = j * 1.0 / (sizes[1]);
				xyz[2] = k * 1.0 / (sizes[2]);
Igor Konshin's avatar
Igor Konshin committed
84
				newverts.push_back(m->CreateNode(xyz)); // Create node in the mesh
85
86
87
				if( ((int)newverts.size()-1) != V_ID(i,j,k))
					printf("v_id = %ld, [%d,%d,%d] = %d\n", newverts.size()-1,i,j,k, V_ID(i, j, k));
			}
Igor Konshin's avatar
Igor Konshin committed
88
		
89
90
91
	std::vector< std::vector<Node *> > c_f_nodes(6);
	for (int i = 0; i < 6; i++)
		c_f_nodes[i].resize(4);
Igor Konshin's avatar
Igor Konshin committed
92

93
94
95
96
97
98
	for(int i = localstart[0]+1; i <= localend[0]; i++)
		for(int j = localstart[1]+1; j <= localend[1]; j++)
			for(int k = localstart[2]+1; k <= localend[2]; k++)
			{
				const INMOST_DATA_ENUM_TYPE nvf[24] = {0,4,6,2,1,3,7,5,0,1,5,4,2,6,7,3,0,2,3,1,4,5,7,6};
				const INMOST_DATA_ENUM_TYPE numnodes[6] = { 4, 4, 4, 4, 4, 4 };
Igor Konshin's avatar
Igor Konshin committed
99
				Node * verts[8];
100
101
102
103
104
105
106
107
				verts[0] = newverts[V_ID(i - 1,j - 1, k - 1)];
				verts[1] = newverts[V_ID(i - 0,j - 1, k - 1)];
				verts[2] = newverts[V_ID(i - 1,j - 0, k - 1)];
				verts[3] = newverts[V_ID(i - 0,j - 0, k - 1)];
				verts[4] = newverts[V_ID(i - 1,j - 1, k - 0)];
				verts[5] = newverts[V_ID(i - 0,j - 1, k - 0)];
				verts[6] = newverts[V_ID(i - 1,j - 0, k - 0)];
				verts[7] = newverts[V_ID(i - 0,j - 0, k - 0)];
Igor Konshin's avatar
Igor Konshin committed
108
109

				m->CreateCell(verts,nvf,numnodes,6).first; // Create the cubic cell in the mesh
110
			}
Igor Konshin's avatar
Igor Konshin committed
111
112
113
114

	m->ResolveShared(); // Resolve duplicate nodes

	m->ExchangeGhost(2,NODE); // Construct Ghost cells in 2 layers connected via nodes
115
116
117
118
119
120
121
122
	return m;
}

Mesh * ParallelCubePrismGenerator(INMOST_MPI_Comm comm, int nx, int ny, int nz)
{
	int procs_per_axis[3] = {1,1,1};
	int sizes[3] = {nx,ny,nz};
	int rank,size;
Igor Konshin's avatar
Igor Konshin committed
123
124
125
	Mesh * m = new Mesh(); // Create a mesh to be constructed

	m->SetCommunicator(comm); // Set the MPI communicator, usually MPI_COMM_WORLD
126
127

#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
128
	//MPI::COMM_WORLD.Set_errhandler(MPI::ERRORS_THROW_EXCEPTIONS);
129
#endif
Igor Konshin's avatar
Igor Konshin committed
130
131
132
133
134

	rank = m->GetProcessorRank(); // Get the rank of the current process
	size = m->GetProcessorsNumber(); // Get the number of processors used in communicator comm

	// Compute the configuration of processes connection
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
	{
		int divsize = size;
		std::vector<int> divs;
		while( divsize > 1 )
		{
			for(int k = 2; k <= divsize; k++)
				if( divsize % k == 0 )
				{
					divs.push_back(k);
					divsize /= k;
					break;
				}
		}
		int elements_per_procs[3] = {sizes[0],sizes[1],sizes[2]};
		for(std::vector<int>::reverse_iterator it = divs.rbegin(); it != divs.rend(); it++)
		{
			int * max = std::max_element(elements_per_procs+0,elements_per_procs+3);
			procs_per_axis[max-elements_per_procs] *= *it;
			(*max) /= *it;
		}
	}
Igor Konshin's avatar
Igor Konshin committed
156

157
	//rank = proc_coords[2] * procs_per_axis[0] *procs_per_axis[1] + proc_coords[1] * procs_per_axis[0] + proc_coords[0];
Igor Konshin's avatar
Igor Konshin committed
158
159
	int proc_coords[3] = {rank % procs_per_axis[0] , rank / procs_per_axis[0] % procs_per_axis[1], rank / (procs_per_axis[0] *procs_per_axis[1]) };

160
	int localsize[3], localstart[3], localend[3];
Igor Konshin's avatar
Igor Konshin committed
161
	int avgsize[3] =
162
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
163
164
165
		(int)ceil((double)sizes[0]/procs_per_axis[0]),
		(int)ceil((double)sizes[1]/procs_per_axis[1]),
		(int)ceil((double)sizes[2]/procs_per_axis[2])
166
	};
Igor Konshin's avatar
Igor Konshin committed
167

168
169
170
171
172
173
174
175
	for(int j = 0; j < 3; j++)
	{
		localstart[j] = avgsize[j] * proc_coords[j];
		if( proc_coords[j] == procs_per_axis[j] - 1 )
			localsize[j] = sizes[j] - avgsize[j] * (procs_per_axis[j]-1);
		else localsize[j] = avgsize[j];
		localend[j] = localstart[j] + localsize[j];
	}
Igor Konshin's avatar
Igor Konshin committed
176

177
178
	std::vector<Node *> newverts;
	newverts.reserve(localsize[0]*localsize[1]*localsize[2]);
Igor Konshin's avatar
Igor Konshin committed
179

180
181
182
183
184
185
186
187
	for(int i = localstart[0]; i <= localend[0]; i++)
		for(int j = localstart[1]; j <= localend[1]; j++)
			for(int k = localstart[2]; k <= localend[2]; k++)
			{
				Storage::real xyz[3];
				xyz[0] = i * 1.0 / (sizes[0]);
				xyz[1] = j * 1.0 / (sizes[1]);
				xyz[2] = k * 1.0 / (sizes[2]);
Igor Konshin's avatar
Igor Konshin committed
188
				newverts.push_back(m->CreateNode(xyz)); // Create node in the mesh
189
190
191
				if( ((int)newverts.size()-1) != V_ID(i,j,k))
					printf("v_id = %ld, [%d,%d,%d] = %d\n", newverts.size()-1,i,j,k, V_ID(i, j, k));
			}
Igor Konshin's avatar
Igor Konshin committed
192

193
194
195
196
197
	std::vector< std::vector<Node *> > c_f_nodes(5);
	for (int i = 0; i < 3; i++)
		c_f_nodes[i].resize(4);
	for (int i = 3; i < 5; i++)
		c_f_nodes[i].resize(3);
Igor Konshin's avatar
Igor Konshin committed
198

199
200
201
202
203
204
	for(int i = localstart[0]+1; i <= localend[0]; i++)
		for(int j = localstart[1]+1; j <= localend[1]; j++)
			for(int k = localstart[2]+1; k <= localend[2]; k++)
			{
				const INMOST_DATA_ENUM_TYPE NE_nvf1[18] = {0,4,6,2,0,3,7,4,2,6,7,3,0,2,3,4,7,6};
				const INMOST_DATA_ENUM_TYPE NE_nvf2[18] = {0,4,7,3,1,3,7,5,0,1,5,4,0,3,1,4,5,7};
Igor Konshin's avatar
Igor Konshin committed
205

206
207
				const INMOST_DATA_ENUM_TYPE NE_nvf3[18] = {0,4,6,2,2,6,5,1,1,5,4,0,0,2,1,4,5,6};
				const INMOST_DATA_ENUM_TYPE NE_nvf4[18] = {1,5,6,2,1,3,7,5,7,3,2,6,1,2,3,6,5,7};
Igor Konshin's avatar
Igor Konshin committed
208

209
				const INMOST_DATA_ENUM_TYPE numnodes[5] = {4,4,4,3,3};
Igor Konshin's avatar
Igor Konshin committed
210
211

				Node * verts[8];
212
213
214
215
216
217
218
219
				verts[0] = newverts[V_ID(i - 1,j - 1, k - 1)];
				verts[1] = newverts[V_ID(i - 0,j - 1, k - 1)];
				verts[2] = newverts[V_ID(i - 1,j - 0, k - 1)];
				verts[3] = newverts[V_ID(i - 0,j - 0, k - 1)];
				verts[4] = newverts[V_ID(i - 1,j - 1, k - 0)];
				verts[5] = newverts[V_ID(i - 0,j - 1, k - 0)];
				verts[6] = newverts[V_ID(i - 1,j - 0, k - 0)];
				verts[7] = newverts[V_ID(i - 0,j - 0, k - 0)];
Igor Konshin's avatar
Igor Konshin committed
220
221

				// Create two prismatic cells in the mesh
222
223
224
225
226
				if ((i + j) % 2 == 0)
				{
					m->CreateCell(verts,NE_nvf1,numnodes,5).first;
					m->CreateCell(verts,NE_nvf2,numnodes,5).first;
				}
Igor Konshin's avatar
Igor Konshin committed
227

228
229
230
231
232
233
				else
				{
					m->CreateCell(verts,NE_nvf3,numnodes,5).first;
					m->CreateCell(verts,NE_nvf4,numnodes,5).first;
				}
			}
Igor Konshin's avatar
Igor Konshin committed
234
235
236
237

	m->ResolveShared(); // Resolve duplicate nodes

	m->ExchangeGhost(2,NODE); // Construct Ghost cells in 2 layers connected via nodes
238
239
240
	return m;
}

Igor Konshin's avatar
Igor Konshin committed
241
int main(int argc, char *argv[])
242
243
244
245
{
	Mesh * mesh;
	Mesh::Initialize(&argc,&argv);

Igor Konshin's avatar
Igor Konshin committed
246
247
248
	// Specify mesh configuration
	int ng = 3, nx = MESH_SIZE, ny = MESH_SIZE, nz = 1, args = 0;
	if( argc > 4 )
249
	{
Igor Konshin's avatar
Igor Konshin committed
250
251
252
253
254
		ng = atoi(argv[1]); // 3 - Prismatic, 4 - Cubic generator
		nx = atoi(argv[2]);
		ny = atoi(argv[3]);
		nz = atoi(argv[4]);
		args = 1;
255
	}
Igor Konshin's avatar
Igor Konshin committed
256
257
258
259
260
261

 	// Construct a mesh in parallel
	//MPI_Barrier(mesh->GetCommunicator());
	double tt = Timer();
	if( ng == 3 )
		mesh = ParallelCubePrismGenerator(INMOST_MPI_COMM_WORLD,nx,ny,nz);
262
	else
Igor Konshin's avatar
Igor Konshin committed
263
		mesh = ParallelCubeGenerator(INMOST_MPI_COMM_WORLD,nx,ny,nz);
264
	tt = Timer() - tt;
Igor Konshin's avatar
Igor Konshin committed
265
	tt = mesh->Integrate(tt)/mesh->GetProcessorsNumber(); //???
266
267
	if( mesh->GetProcessorRank() == 0 )
	{
Igor Konshin's avatar
Igor Konshin committed
268
269
270
271
272
273
		if( args == 0 ) std::cout << "Default ";
		if( ng == 3 )   std::cout << "Prismatic ";
		else            std::cout << "Cubic ";
		std::cout << "Grid: " << nx << " x " << ny << " x " << nz << std::endl;
		std::cout << "Processors: " <<mesh->GetProcessorsNumber() << std::endl;	
		std::cout << "Mesh generator time: " << tt << std::endl;
274
	}
Igor Konshin's avatar
Igor Konshin committed
275
276
277
278
279
280
281
282
283
284
285
286
287

	std::string filename = "grid";
	if( mesh->GetProcessorsNumber() == 1 )
		filename += ".vtk";
	else
		filename += ".pvtk";
	MPI_Barrier(mesh->GetCommunicator());
	tt = Timer();
	mesh->Save(filename); // Save constructed mesh to the file
	MPI_Barrier(mesh->GetCommunicator());
	tt = Timer() - tt;
	if( mesh->GetProcessorRank() == 0 ) std::cout << "Save to file \"" << filename << "\" time: " << tt << std::endl;

288
289
290
291
	delete mesh;
	Mesh::Finalize();
	return 0;
}