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

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

9
#define MESH_SIZE 32
10
11
12

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

13
14
15
16
17
18
19
20
21
22
23
24
25
/*      (4)*-------*(6)
          /|      /|
         /       / |
        /  |    /  |
    (5)*-------*(7)|
       |   |   |   |
       |       |   |
       |   |   |   |
       |(0)*- -|- -*(2)
       |  /    |  /
       |       | /
       |/      |/
    (1)*-------*(3)      */
26
27
void CreateCubeElement(Mesh *m, ElementArray<Node> verts)
{
28
29
30
	// Define six cube faces assuming verts are numerated in the way presented above
	const INMOST_DATA_INTEGER_TYPE face_nodes[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_INTEGER_TYPE num_nodes[6]   = {4,       4,       4,       4,       4,       4};
Igor Konshin's avatar
Igor Konshin committed
31

32
33
	m->CreateCell(verts,face_nodes,num_nodes,6); // Create the cubic cell in the mesh
}
Igor Konshin's avatar
Igor Konshin committed
34

35
36
37
38
39
40
41
42
43
44
45
46
47
/*      (4)*-------*(6)
          /|  __/ /|
         / __/   / |
        / /|    /  |
    (5)*-------*(7)|
       |   |   |   |
       |       |   |
       |   |   |   |
       |(0)*- -|- -*(2)
       |  /    |/ /
       |    _/ | /
       |/_/    |/
    (1)*-------*(3)      */
48
49
void CreateNEPrismElements(Mesh *m, ElementArray<Node> verts)
{
50
51
52
53
54
55
56
57
58
	// Define prism faces assuming verts are numerated in the way presented above
	const INMOST_DATA_INTEGER_TYPE nw_face_nodes[18] = {0,4,6,2, 2,6,5,1, 1,5,4,0, 0,2,1, 4,5,6};
	const INMOST_DATA_INTEGER_TYPE nw_num_nodes[5]   = {4,       4,       4,       3,     3};
	const INMOST_DATA_INTEGER_TYPE se_face_nodes[18] = {1,5,6,2, 1,3,7,5, 7,3,2,6, 1,2,3, 6,5,7};
	const INMOST_DATA_INTEGER_TYPE se_num_nodes[5]   = {4,       4,       4,       3,     3};

	m->CreateCell(verts,nw_face_nodes,nw_num_nodes,5); // Create north-west prismatic cell
	m->CreateCell(verts,se_face_nodes,se_num_nodes,5); // Create south-east prismatic cell
}
Igor Konshin's avatar
Igor Konshin committed
59

60
61
62
63
64
65
66
67
68
69
70
71
72
/*      (4)*-------*(6)
          /|\     /|
         /   \   / |
        /  |  \ /  |
    (5)*-------*(7)|
       |   |   |   |
       |       |   |
       |   |   |   |
       |(0)*- -|- -*(2)
       |  / \  |  /
       |       | /
       |/     \|/
    (1)*-------*(3)      */
73
74
void CreateNWPrismElements(Mesh *m, ElementArray<Node> verts)
{
75
76
77
78
79
80
81
82
	// Define prism faces assuming verts are numerated in the way presented above
	const INMOST_DATA_INTEGER_TYPE ne_face_nodes[18] = {0,4,6,2, 0,3,7,4, 2,6,7,3, 0,2,3, 4,7,6};
	const INMOST_DATA_INTEGER_TYPE ne_num_nodes[5]   = {4,       4,       4,       3,     3};
	const INMOST_DATA_INTEGER_TYPE sw_face_nodes[18] = {0,4,7,3, 1,3,7,5, 0,1,5,4, 0,3,1, 4,5,7};
	const INMOST_DATA_INTEGER_TYPE sw_num_nodes[5]   = {4,       4,       4,       3,     3};

	m->CreateCell(verts,ne_face_nodes,ne_num_nodes,5); // Create north-east prismatic cell
	m->CreateCell(verts,sw_face_nodes,sw_num_nodes,5); // Create south-west prismatic cell
Dmitry Bagaev's avatar
Dmitry Bagaev committed
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
}


/*      (4)*-------*(6)
          /|\     /|
         /   \   / |
        /  |  \ /  |
    (5)*-------*(7)|
       |   |   |   |
       |       |   |
       |   |   |   |
       |(0)*- -|- -*(2)
       |  / \  |  /
       |       | /
       |/     \|/
    (1)*-------*(3)      */
99
100
void CreateNWTetElements(Mesh *m, ElementArray<Node> verts)
{
ahmadabushaikha's avatar
ahmadabushaikha committed
101
102
103
	// Define prism faces assuming verts are numerated in the way presented above
	const INMOST_DATA_INTEGER_TYPE ne_face_nodes1[12] = {0,1,5,  5,1,3,   1,0,3, 3,0,5};
	const INMOST_DATA_INTEGER_TYPE ne_num_nodes1[4]   = {3,3,3,3};
Dmitry Bagaev's avatar
Dmitry Bagaev committed
104

ahmadabushaikha's avatar
ahmadabushaikha committed
105
106
	const INMOST_DATA_INTEGER_TYPE ne_face_nodes2[12] = {0,3,5, 0,7,3, 5,3,7, 0,5,7};
	const INMOST_DATA_INTEGER_TYPE ne_num_nodes2[4]   = {3,3,3,3};
Dmitry Bagaev's avatar
Dmitry Bagaev committed
107

ahmadabushaikha's avatar
ahmadabushaikha committed
108
109
	const INMOST_DATA_INTEGER_TYPE ne_face_nodes3[12] = {0,7,5, 4,5,7, 0,5,4, 0,4,7};
	const INMOST_DATA_INTEGER_TYPE ne_num_nodes3[4]   = {3,3,3,3};
Dmitry Bagaev's avatar
Dmitry Bagaev committed
110
111


ahmadabushaikha's avatar
ahmadabushaikha committed
112
113
	const INMOST_DATA_INTEGER_TYPE sw_face_nodes1[12] = {0,3,7, 2,7,3, 0,7,2, 0,2,3};
	const INMOST_DATA_INTEGER_TYPE sw_num_nodes1[4]   = {3,3,3,3};
Dmitry Bagaev's avatar
Dmitry Bagaev committed
114

ahmadabushaikha's avatar
ahmadabushaikha committed
115
116
	const INMOST_DATA_INTEGER_TYPE sw_face_nodes2[12] = {0,7,4, 0,2,7, 2,4,7, 0,4,2};
	const INMOST_DATA_INTEGER_TYPE sw_num_nodes2[4]   = {3,3,3,3};
Dmitry Bagaev's avatar
Dmitry Bagaev committed
117

ahmadabushaikha's avatar
ahmadabushaikha committed
118
119
	const INMOST_DATA_INTEGER_TYPE sw_face_nodes3[12] = {4,6,2, 6,7,2, 4,7,6, 4,2,7};
	const INMOST_DATA_INTEGER_TYPE sw_num_nodes3[4]   = {3,3,3,3};
Dmitry Bagaev's avatar
Dmitry Bagaev committed
120

ahmadabushaikha's avatar
ahmadabushaikha committed
121
122
123
124
125
126
	m->CreateCell(verts,ne_face_nodes1,ne_num_nodes1,4); // Create north-east prismatic cell
	m->CreateCell(verts,ne_face_nodes2,ne_num_nodes2,4); // Create north-east prismatic cell
	m->CreateCell(verts,ne_face_nodes3,ne_num_nodes3,4); // Create north-east prismatic cell
	m->CreateCell(verts,sw_face_nodes1,sw_num_nodes1,4); // Create south-west prismatic cell
	m->CreateCell(verts,sw_face_nodes2,sw_num_nodes2,4); // Create south-west prismatic cell
	m->CreateCell(verts,sw_face_nodes3,sw_num_nodes3,4); // Create south-west prismatic cell
127
128
}

129
130
Mesh * ParallelGenerator(INMOST_MPI_Comm comm, int ng, int nx, int ny, int nz)
{
131
132
133
	int procs_per_axis[3] = {1,1,1};
	int sizes[3] = {nx,ny,nz};
	int rank,size;
Igor Konshin's avatar
Igor Konshin committed
134
135
136
	Mesh * m = new Mesh(); // Create a mesh to be constructed

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

#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
139
	MPI_Comm_set_errhandler(comm,MPI_ERRORS_RETURN);
Kirill Terekhov's avatar
Kirill Terekhov committed
140
	//MPI::COMM_WORLD.Set_errhandler(MPI::ERRORS_THROW_EXCEPTIONS);
141
#endif
Igor Konshin's avatar
Igor Konshin committed
142
143
144
145
146

	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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
	{
		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
168

169
	//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
170
171
	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]) };

172
	int localsize[3], localstart[3], localend[3];
Igor Konshin's avatar
Igor Konshin committed
173
	int avgsize[3] =
174
175
176
177
178
			{
					(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])
			};
Igor Konshin's avatar
Igor Konshin committed
179

180
181
182
183
184
185
186
187
	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
188

189
	// Create i-j-k structure of nodes
190
	ElementArray<Node> newverts(m);
191
	newverts.reserve(localsize[0]*localsize[1]*localsize[2]);
Igor Konshin's avatar
Igor Konshin committed
192

193
194
195
196
197
198
199
200
	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]);
201
				newverts.push_back(m->CreateNode(xyz)); // Create node in the mesh with index V_ID(i,j,k)
202
			}
Igor Konshin's avatar
Igor Konshin committed
203

204
	// Create i-j-k structure of elements
205
206
207
208
	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++)
			{
209
210
211
212
213
214
215
216
217
218
219
				// Create local array of eight nodes                           /*      (4)*-------*(6)  */
				// using representation on the right figure                    /*        /|      /|     */
				ElementArray<Node> verts(m);                                   /*       /       / |     */
				verts.push_back(newverts[V_ID(i - 1,j - 1, k - 1)]); // 0      /*      /  |    /  |     */
				verts.push_back(newverts[V_ID(i - 0,j - 1, k - 1)]); // 1      /*  (5)*-------*(7)|     */
				verts.push_back(newverts[V_ID(i - 1,j - 0, k - 1)]); // 2      /*     |   |   |   |     */
				verts.push_back(newverts[V_ID(i - 0,j - 0, k - 1)]); // 3      /*     |       |   |     */
				verts.push_back(newverts[V_ID(i - 1,j - 1, k - 0)]); // 4      /*     |   |   |   |     */
				verts.push_back(newverts[V_ID(i - 0,j - 1, k - 0)]); // 5      /*     |(0)*- -|- -*(2)  */
				verts.push_back(newverts[V_ID(i - 1,j - 0, k - 0)]); // 6      /*     |  /    |  /      */
				verts.push_back(newverts[V_ID(i - 0,j - 0, k - 0)]); // 7      /*     |       | /       */
220
				/*     |/      |/        */
221
				// Create cells based on parameter ng                          /*  (1)*-------*(3)      */
Kirill Terekhov's avatar
Kirill Terekhov committed
222
				if (ng == 5) // Create tetrahedral grid
ahmadabushaikha's avatar
ahmadabushaikha committed
223
224
225
226
				{
					CreateNWTetElements(m,verts);
				}
				else if (ng == 4) // Create cubic cell
227
				{
228
					CreateCubeElement(m,verts);
229
				}
230
231
232
233
234
				else if ((i + j) % 2 == 0) // Create two prism cells
				{
					CreateNWPrismElements(m,verts);
				}
				else // Two prism cells with different diagonal direction
235
				{
236
					CreateNEPrismElements(m,verts);
237
238
				}
			}
Igor Konshin's avatar
Igor Konshin committed
239
240
241
242

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

	m->ExchangeGhost(2,NODE); // Construct Ghost cells in 2 layers connected via nodes
243
	return m;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
244
245
}

246
247
int main(int argc, char *argv[])
{
248
	std::string mesh_name;
249
250
251
	Mesh * mesh;
	Mesh::Initialize(&argc,&argv);

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

263
	// Construct a mesh in parallel
Igor Konshin's avatar
Igor Konshin committed
264
265
	//MPI_Barrier(mesh->GetCommunicator());
	double tt = Timer();
266
	mesh = ParallelGenerator(INMOST_MPI_COMM_WORLD,ng,nx,ny,nz);
267
	tt = Timer() - tt;
Igor Konshin's avatar
Igor Konshin committed
268
	tt = mesh->Integrate(tt)/mesh->GetProcessorsNumber(); //???
269
270
	if( mesh->GetProcessorRank() == 0 )
	{
271
		if( args == 0 ) std::cout << "Usage: " << argv[0] << " ng nx ny nz [output.[p]vtk]" << std::endl << "ng - 3 for prismatic mesh, 4 for cubic mesh" << std::endl;
Igor Konshin's avatar
Igor Konshin committed
272
		if( args == 0 ) std::cout << "Default ";
Kirill Terekhov's avatar
Kirill Terekhov committed
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
		if( ng == 5 )
		{
			std::cout << "Tetrahedral ";
			std::stringstream str;
			str << "TETRAHEDRAL_" << nx << "x" << ny << "x" << nz;
			mesh_name = str.str();
		}
		else if( ng == 4 )
		{
			std::cout << "Cubic ";
			std::stringstream str;
			str << "CUBIC_" << nx << "x" << ny << "x" << nz;
			mesh_name = str.str();
		}
		else
		{
			std::cout << "Prismatic ";
			std::stringstream str;
			str << "PRISMATIC_" << nx << "x" << ny << "x" << nz;
			mesh_name = str.str();
		}
Igor Konshin's avatar
Igor Konshin committed
294
		std::cout << "Grid: " << nx << " x " << ny << " x " << nz << std::endl;
295
		std::cout << "Processors: " <<mesh->GetProcessorsNumber() << std::endl;
Igor Konshin's avatar
Igor Konshin committed
296
		std::cout << "Mesh generator time: " << tt << std::endl;
297
	}
Igor Konshin's avatar
Igor Konshin committed
298

299
300
301
302
303
	std::string filename;
	if (argc > 5)
	{
		filename = argv[5];
	}
Igor Konshin's avatar
Igor Konshin committed
304
	else
305
306
307
308
309
310
311
	{
		filename = "grid";
		if( mesh->GetProcessorsNumber() == 1 )
			filename += ".vtk";
		else
			filename += ".pvtk";
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
312

313
314
	Storage::bulk_array name = mesh->self()->BulkArray(mesh->CreateTag("GRIDNAME",DATA_BULK,MESH,NONE));
	name.replace(name.begin(),name.end(),mesh_name.begin(),mesh_name.end());
Kirill Terekhov's avatar
Kirill Terekhov committed
315

316
#if defined(USE_MPI)
Igor Konshin's avatar
Igor Konshin committed
317
	MPI_Barrier(mesh->GetCommunicator());
318
#endif
Igor Konshin's avatar
Igor Konshin committed
319
320
	tt = Timer();
	mesh->Save(filename); // Save constructed mesh to the file
321
#if defined(USE_MPI)
Igor Konshin's avatar
Igor Konshin committed
322
	MPI_Barrier(mesh->GetCommunicator());
323
#endif
Igor Konshin's avatar
Igor Konshin committed
324
325
326
	tt = Timer() - tt;
	if( mesh->GetProcessorRank() == 0 ) std::cout << "Save to file \"" << filename << "\" time: " << tt << std::endl;

327
328
329
	delete mesh;
	Mesh::Finalize();
	return 0;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
330
}