main.cpp 13.4 KB
Newer Older
1
2
3
4
#include "amesh.h"

using namespace INMOST;

Kirill Terekhov's avatar
Kirill Terekhov committed
5
bool output_file = true;
Kirill Terekhov's avatar
Kirill Terekhov committed
6
7
8
bool balance_mesh = true;
bool balance_mesh_refine = false;
bool balance_mesh_coarse = false;
9
std::string file_format = ".pvtu";
10

11
#ifndef M_PI
Kirill Terekhov's avatar
Kirill Terekhov committed
12
#define M_PI 3.1415926535897932
13
14
#endif

15
16
17
18
19
20
int main(int argc, char ** argv)
{
	Mesh::Initialize(&argc,&argv);
	
	if( argc > 1 )
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
21
		Mesh m;
22
		
23
		m.SetCommunicator(INMOST_MPI_COMM_WORLD);
Kirill Terekhov's avatar
Kirill Terekhov committed
24
25
26
27
28
		if( m.isParallelFileFormat(argv[1]) )
			m.Load(argv[1]);
		else if( m.GetProcessorRank() == 0 )
			m.Load(argv[1]);
	
29
30
31
32
33
34
35
36
		//~ m.SetTopologyCheck(PROHIBIT_MULTIPOLYGON | PROHIBIT_MULTILINE 
						   //~ | PROHIBIT_CONCAVE_CELL
						   //~ | FACE_EDGES_ORDER
						   //~ | DEGENERATE_CELL | DEGENERATE_EDGE | DEGENERATE_FACE
						   //~ | PRINT_NOTIFY | TRIPLE_SHARED_FACE | FLATTENED_CELL
						   //~ | INTERLEAVED_FACES | NEED_TEST_CLOSURE
						   //~ | DUPLICATE_EDGE | DUPLICATE_FACE | DUPLICATE_CELL
						   //~ | ADJACENT_DUPLICATE | ADJACENT_HIDDEN | ADJACENT_VALID | ADJACENT_DIMENSION);
37
		//m.RemTopologyCheck(THROW_EXCEPTION);
38
39
40
41
42
43
44
45
		
		{
			Mesh::GeomParam table;
			table[ORIENTATION] = FACE;
			table[NORMAL] = FACE;
			table[CENTROID] = CELL | FACE;
			m.PrepareGeometricData(table);
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
46
#if defined(USE_PARTITIONER)
47
		std::vector<Storage::integer> nc(m.GetProcessorsNumber());
48
		Partitioner p(&m);
Kirill Terekhov's avatar
Kirill Terekhov committed
49
50
		if( true )
		{
51
			std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "init before "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
52
53
			//p.SetMethod(Partitioner::INNER_KMEANS,Partitioner::Partition);
			p.SetMethod(Partitioner::Parmetis,Partitioner::Partition);
54
			//p.SetMethod(Partitioner::Parmetis,Partitioner::Refine);
Kirill Terekhov's avatar
Kirill Terekhov committed
55
56
57
			p.Evaluate();
			m.Redistribute();
			m.ReorderEmpty(CELL|FACE|EDGE|NODE);
58
59
60
			
			std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "init after "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
			m.Barrier();
61
			p.SetMethod(Partitioner::Parmetis,Partitioner::Repartition);
Kirill Terekhov's avatar
Kirill Terekhov committed
62
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
63
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
64
		m.ExchangeGhost(1,FACE);
Kirill Terekhov's avatar
Kirill Terekhov committed
65
		AdaptiveMesh am(m);
66
67
68
69
		//m.SetTopologyCheck(NEED_TEST_CLOSURE);
		//m.SetTopologyCheck(PROHIBIT_MULTILINE);
		//m.SetTopologyCheck(PROHIBIT_MULTIPOLYGON);
		TagInteger indicator = m.CreateTag("INDICATOR",DATA_INTEGER,CELL,NONE,1);
70
71
72
73
		TagReal wgt = m.CreateTag("WEIGHT",DATA_REAL,CELL,NONE,1);
#if defined(USE_PARTITIONER)
		p.SetWeight(wgt);
#endif
74
75
76

		//m.Save("init"+file_format);
		
Kirill Terekhov's avatar
Kirill Terekhov committed
77
78
79
80
81
82
83
84
		/*
		for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
			indicator[*it] = 1;
		if( !m.Refine(indicator) ) std::cout << "refine failed" << std::endl;
		else std::cout << "refine ok!" << std::endl;
		m.Save("grid.pmf");
		return 0;
		*/
85
86
87
88
89
		int max_levels = 2;
		if( argc > 2 ) max_levels = atoi(argv[2]);
		//bounding box around mesh
		Storage::real cmax[3] = {-1.0e20,-1.0e20,-1.0e20};
		Storage::real cmin[3] = {+1.0e20,+1.0e20,+1.0e20};
90
		Storage::real xyz[3], r, q, cnt[3], cnt0[3], r0;
91
92
93
94
95
96
97
98
99
		//find bounding box around mesh
		for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it)
		{
			for(int d = 0; d < 3; ++d)
			{
				if( it->Coords()[d] > cmax[d] ) cmax[d] = it->Coords()[d];
				if( it->Coords()[d] < cmin[d] ) cmin[d] = it->Coords()[d];
			}
		}
100
101
		m.AggregateMax(cmax,3);
		m.AggregateMin(cmin,3);
102
		r0 = 1;
103
104
		for(int d = 0; d < 3; ++d)
		{
105
			r0 *= cmax[d]-cmin[d];
106
			cnt[d] = (cmax[d]+cmin[d])*0.5;
107
			cnt0[d] = cnt[d];
108
		}
109
110
111
		//r = pow(r,1.0/3.0)/20.0;
		r0 = pow(r0,1.0/3.0);
		r = r0/8.0;
112
		
113
		//int ncells, nfaces, nedges, nnodes;
Kirill Terekhov's avatar
sync    
Kirill Terekhov committed
114
115

		double time = Timer();
116
		
Kirill Terekhov's avatar
Kirill Terekhov committed
117
		for(int k = 0; k < 80; ++k)
118
		{
119

Kirill Terekhov's avatar
sync    
Kirill Terekhov committed
120

Kirill Terekhov's avatar
Kirill Terekhov committed
121
			double step_time = Timer(), ref_time, crs_time, redist_time = 0, tmp_time;
Kirill Terekhov's avatar
sync    
Kirill Terekhov committed
122

Kirill Terekhov's avatar
Kirill Terekhov committed
123
124
			cnt[0] = cnt0[0] + 0.25*r0*sin(k/20.0*M_PI);
			cnt[1] = cnt0[1] + 0.25*r0*cos(k/20.0*M_PI);
125

126
			m.ClearFile();
127
			
128
			std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "start "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
129
			
Kirill Terekhov's avatar
Kirill Terekhov committed
130
			ref_time = Timer();
131
			Storage::integer numref;
132
133
134
135
			int refcnt = 0;
			do
			{
				numref = 0;
136
				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) 
137
				{
138
					if( am.GetLevel(it->self()) < max_levels )
139
					{
140
141
142
						it->Barycenter(xyz);
						//refine a circle
						q = 0;
143
						for(int d = 0; d < 2; ++d)
144
145
							q += (xyz[d]-cnt[d])*(xyz[d]-cnt[d]);
						q = sqrt(q);
146
147
						//if( q < r*(k+1) && q > r*k)
						if( q < r*1.4 && q > r)
148
149
150
151
						{
							indicator[*it] = 1;
							numref++;
						}
152
					}
153
					else indicator[*it] = 0;
154
				}
155
156
				numref = m.Integrate(numref);
				if( numref )
157
				{
158
					//~ std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "refine_" << k << "_" << refcnt << " "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
Kirill Terekhov's avatar
Kirill Terekhov committed
159
					tmp_time = Timer();
160
161
162
#if defined(USE_PARTITIONER)
					if( balance_mesh_refine && refcnt == 0)
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
163
						//m.Barrier();
164
165
						am.ComputeWeightRefine(indicator,wgt);
						p.SetWeight(wgt);
Kirill Terekhov's avatar
Kirill Terekhov committed
166
						//std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "refine before "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
Kirill Terekhov's avatar
Kirill Terekhov committed
167
						if( m.GetProcessorRank() == 0 ) std::cout << "call Evaluate refine" << std::endl;
168
						p.Evaluate();
Kirill Terekhov's avatar
Kirill Terekhov committed
169
						if( m.GetProcessorRank() == 0 ) std::cout << "call Redistribute refine" << std::endl;
170
						m.Redistribute();
Kirill Terekhov's avatar
Kirill Terekhov committed
171
						if( m.GetProcessorRank() == 0 ) std::cout << "balance done refine" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
172
						//std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "refine after "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
Kirill Terekhov's avatar
Kirill Terekhov committed
173
						//m.Barrier();
174
175
					}
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
176
177
178
179
180
181
182
					tmp_time = Timer() - tmp_time;
					redist_time += tmp_time;
					ref_time += tmp_time;
					//ncells = m.TotalNumberOf(CELL);
					//nfaces = m.TotalNumberOf(FACE);
					//nedges = m.TotalNumberOf(EDGE);
					//nnodes = m.TotalNumberOf(NODE);
183
					if( m.GetProcessorRank() == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
184
						std::cout << "beg k " << k << " refcnt " << refcnt << std::endl;// << " cells " << ncells << " faces " << nfaces << " nedges " << nedges << " nnodes " << nnodes << std::endl;
185
186
187
188
					//m.BeginSequentialCode();
					//std::cout << m.GetProcessorRank() << " cells " << m.NumberOfCells() << std::endl;
					//m.EndSequentialCode();

Kirill Terekhov's avatar
Kirill Terekhov committed
189
					if (!am.Refine(indicator)) break;
Kirill Terekhov's avatar
Kirill Terekhov committed
190
					
191
					
Kirill Terekhov's avatar
Kirill Terekhov committed
192
193
194
195
					//ncells = m.TotalNumberOf(CELL);
					//nfaces = m.TotalNumberOf(FACE);
					//nedges = m.TotalNumberOf(EDGE);
					//nnodes = m.TotalNumberOf(NODE);
196
					if( m.GetProcessorRank() == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
197
						std::cout << "end k " << k << " refcnt " << refcnt << std::endl;// << " cells " << ncells << " faces " << nfaces << " nedges " << nedges << " nnodes " << nnodes << std::endl;
198
					
Kirill Terekhov's avatar
Kirill Terekhov committed
199
200
					if( false )
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
201
						std::stringstream file;
202
						file << "ref_" << k << "_" << refcnt << file_format;
Kirill Terekhov's avatar
Kirill Terekhov committed
203
204
						m.Save(file.str());
						if( m.GetProcessorRank() == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
205
							std::cout << "Save " << file.str() << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
206
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
207
					
208
209
210
211
				}
				refcnt++;
			}
			while(numref);
Kirill Terekhov's avatar
Kirill Terekhov committed
212
213
214
			ref_time = Timer() - ref_time;
			
			crs_time = Timer();
215
216
217
218
			refcnt = 0;
			do
			{
				numref = 0;
219
				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
220
				{
221
					if( am.GetLevel(it->self()) > 0 )
222
					{
223
224
225
						it->Barycenter(xyz);
						//refine a circle
						q = 0;
226
						for(int d = 0; d < 2; ++d)
227
228
							q += (xyz[d]-cnt[d])*(xyz[d]-cnt[d]);
						q = sqrt(q);
229
230
						//if( q < r*k)
						if( !(q < r*1.4 && q > r) )
231
232
233
234
						{
							indicator[*it] = 1;
							numref++;
						}
235
					}
236
					else indicator[*it] = 0;
237
				}
238
239
				numref = m.Integrate(numref);
				if( numref )
240
				{
241
					//~ std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "coarse_" << k << "_" << refcnt << " "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
Kirill Terekhov's avatar
Kirill Terekhov committed
242
					tmp_time = Timer();
243
244
245
246
247
248
#if defined(USE_PARTITIONER)
					if( balance_mesh_coarse && refcnt == 0)
					{
						m.Barrier();
						am.ComputeWeightCoarse(indicator,wgt);
						p.SetWeight(wgt);
Kirill Terekhov's avatar
Kirill Terekhov committed
249
						//std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "coarse before "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
Kirill Terekhov's avatar
Kirill Terekhov committed
250
						if( m.GetProcessorRank() == 0 ) std::cout << "call Evaluate coarse" << std::endl;
251
						p.Evaluate();
Kirill Terekhov's avatar
Kirill Terekhov committed
252
						if( m.GetProcessorRank() == 0 ) std::cout << "call Redistribute coarse" << std::endl;
253
						m.Redistribute();
Kirill Terekhov's avatar
Kirill Terekhov committed
254
						if( m.GetProcessorRank() == 0 ) std::cout << "balance done coarse" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
255
						//std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "coarse after "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
256
257
258
						m.Barrier();
					}
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
259
260
261
					tmp_time = Timer() - tmp_time;
					redist_time += tmp_time;
					crs_time += tmp_time;
262
					
Kirill Terekhov's avatar
Kirill Terekhov committed
263
264
265
266
					//ncells = m.TotalNumberOf(CELL);
					//nfaces = m.TotalNumberOf(FACE);
					//nedges = m.TotalNumberOf(EDGE);
					//nnodes = m.TotalNumberOf(NODE);
267
					if( m.GetProcessorRank() == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
268
						std::cout << ":beg k " << k << " crscnt " << refcnt << std::endl;//" cells " << ncells << " faces " << nfaces << " nedges " << nedges << " nnodes " << nnodes << std::endl;
269
270
271
					//m.BeginSequentialCode();
					//std::cout << m.GetProcessorRank() << " cells " << m.NumberOfCells() << std::endl;
					//m.EndSequentialCode();
Kirill Terekhov's avatar
Kirill Terekhov committed
272
273
					
					if( !am.Coarse(indicator) ) break;
Kirill Terekhov's avatar
Kirill Terekhov committed
274
					
Kirill Terekhov's avatar
Kirill Terekhov committed
275
276
277
278
					//ncells = m.TotalNumberOf(CELL);
					//nfaces = m.TotalNumberOf(FACE);
					//nedges = m.TotalNumberOf(EDGE);
					//nnodes = m.TotalNumberOf(NODE);
279
					if( m.GetProcessorRank() == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
280
						std::cout << ":end k " << k << " crscnt " << refcnt << std::endl;//" cells " << ncells << " faces " << nfaces << " nedges " << nedges << " nnodes " << nnodes << std::endl;
281
					
Kirill Terekhov's avatar
Kirill Terekhov committed
282
283
284
					if( false ) 
					{
						std::stringstream file;
285
						file << "crs_" << k << "_" << refcnt << file_format;
Kirill Terekhov's avatar
Kirill Terekhov committed
286
287
288
289
290
						m.Save(file.str());
						if( m.GetProcessorRank() == 0 )
							std::cout << "Save " << file.str() << std::endl;
					}
					
291
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
292
				
293
294
295
				refcnt++;
			}
			while(numref);
Kirill Terekhov's avatar
Kirill Terekhov committed
296
			crs_time = Timer() - crs_time;
297
			
Kirill Terekhov's avatar
Kirill Terekhov committed
298
			tmp_time = Timer();
299
#if defined(USE_PARTITIONER)
300
			if( balance_mesh )
301
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
302
				//m.Barrier();
303
				p.SetWeight(Tag());
304
				//~ std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "finish before "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
Kirill Terekhov's avatar
Kirill Terekhov committed
305
				if( m.GetProcessorRank() == 0 ) std::cout << "call Evaluate" << std::endl;
306
				p.Evaluate();
Kirill Terekhov's avatar
Kirill Terekhov committed
307
				if( m.GetProcessorRank() == 0 ) std::cout << "call Redistribute" << std::endl;
308
				m.Redistribute();
Kirill Terekhov's avatar
Kirill Terekhov committed
309
				if( m.GetProcessorRank() == 0 ) std::cout << "balance done" << std::endl;
310
				//~ std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "finish after "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
Kirill Terekhov's avatar
Kirill Terekhov committed
311
				//m.Barrier();
312
313
			}
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
314
			tmp_time = Timer() - tmp_time;
Kirill Terekhov's avatar
Kirill Terekhov committed
315
			redist_time += tmp_time;
316
			
317
			std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "finish "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
Kirill Terekhov's avatar
Kirill Terekhov committed
318
			 
319
320
321
			 

			
322
			if( output_file )
323
			{
324
325
326
327
328
329
330
331
				TagInteger tag_owner = m.CreateTag("OWN",DATA_INTEGER,CELL,NONE,1);
				TagInteger tag_owner0 = m.GetTag("OWNER_PROCESSOR");
				TagInteger tag_stat = m.CreateTag("STAT",DATA_INTEGER,CELL,NONE,1);
				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
				{
					tag_owner[*it] = tag_owner0[*it];
					tag_stat[*it] = it->GetStatus();
				}
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
				{
					std::stringstream file;
					file << "step_" << k << file_format;
					m.Save(file.str());
					if( m.GetProcessorRank() == 0 )
						std::cout << "Save " << file.str() << std::endl;
				}
				if( false )
				{
					std::stringstream file;
					file << "step_" << k << ".pmf";
					m.Save(file.str());
					if( m.GetProcessorRank() == 0 )
						std::cout << "Save " << file.str() << std::endl;
				}
347
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
348
			else if( m.GetProcessorRank() == 0 ) std::cout << "step " << k << std::endl;
Kirill Terekhov's avatar
sync    
Kirill Terekhov committed
349
350
			
			step_time = Timer() - step_time;
Kirill Terekhov's avatar
Kirill Terekhov committed
351
			if( m.GetProcessorRank() == 0 ) std::cout << "step time " << step_time << " refine " << ref_time << " coarse " << crs_time << " balance " << redist_time << std::endl;
352
353
			//~ MPI_Barrier(MPI_COMM_WORLD);
			//~ MPI_Abort(MPI_COMM_WORLD,-1);
354
		}
Kirill Terekhov's avatar
sync    
Kirill Terekhov committed
355
356
357
358


		time = Timer() - time;
		if( m.GetProcessorRank() == 0 ) std::cout << "total time: " << time << " processors " << m.GetProcessorsNumber() << std::endl;
359
360
	}
	else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2]" << std::endl;
361
362
	
	Mesh::Finalize();
Kirill Terekhov's avatar
Kirill Terekhov committed
363
	return 0;
364
}