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

using namespace INMOST;

5
6
std::string file_format = ".pmf";

7
8
9
10
11
12
int main(int argc, char ** argv)
{
	Mesh::Initialize(&argc,&argv);
	
	if( argc > 1 )
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
13
		Mesh m;
14
		
15
		m.SetCommunicator(INMOST_MPI_COMM_WORLD);
Kirill Terekhov's avatar
Kirill Terekhov committed
16
17
18
19
20
		if( m.isParallelFileFormat(argv[1]) )
			m.Load(argv[1]);
		else if( m.GetProcessorRank() == 0 )
			m.Load(argv[1]);
	
21
22
		//m.SetTopologyCheck(PROHIBIT_MULTIPOLYGON | PROHIBIT_MULTILINE | DEGENERATE_CELL | DEGENERATE_EDGE | DEGENERATE_FACE | PRINT_NOTIFY | TRIPLE_SHARED_FACE | FLATTENED_CELL | INTERLEAVED_FACES | NEED_TEST_CLOSURE);
		//m.RemTopologyCheck(THROW_EXCEPTION);
Kirill Terekhov's avatar
Kirill Terekhov committed
23
#if defined(USE_PARTITIONER)
24
		Partitioner p(&m);
Kirill Terekhov's avatar
Kirill Terekhov committed
25
26
		if( true )
		{
27
			m.Barrier();
28
			std::cout << "before on " << m.GetProcessorRank() << " " << m.NumberOfCells() << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
29
			p.SetMethod(Partitioner::INNER_KMEANS,Partitioner::Partition);
30
			//p.SetMethod(Partitioner::Parmetis,Partitioner::Refine);
Kirill Terekhov's avatar
Kirill Terekhov committed
31
32
33
			p.Evaluate();
			m.Redistribute();
			m.ReorderEmpty(CELL|FACE|EDGE|NODE);
34
			std::cout << "after on " << m.GetProcessorRank() << " " << m.NumberOfCells() << std::endl;
35
			p.SetMethod(Partitioner::Parmetis,Partitioner::Repartition);
Kirill Terekhov's avatar
Kirill Terekhov committed
36
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
37
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
38
		m.ExchangeGhost(2,FACE);
Kirill Terekhov's avatar
Kirill Terekhov committed
39
		AdaptiveMesh am(m);
40
41
42
43
		//m.SetTopologyCheck(NEED_TEST_CLOSURE);
		//m.SetTopologyCheck(PROHIBIT_MULTILINE);
		//m.SetTopologyCheck(PROHIBIT_MULTIPOLYGON);
		TagInteger indicator = m.CreateTag("INDICATOR",DATA_INTEGER,CELL,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
44
45
46
47
48
49
50
51
		/*
		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;
		*/
52
53
54
55
56
		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};
57
		Storage::real xyz[3], r, q, cnt[3], cnt0[3], r0;
58
59
60
61
62
63
64
65
66
		//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];
			}
		}
67
68
		m.AggregateMax(cmax,3);
		m.AggregateMin(cmin,3);
69
		r0 = 1;
70
71
		for(int d = 0; d < 3; ++d)
		{
72
			r0 *= cmax[d]-cmin[d];
73
			cnt[d] = (cmax[d]+cmin[d])*0.5;
74
			cnt0[d] = cnt[d];
75
		}
76
77
78
		//r = pow(r,1.0/3.0)/20.0;
		r0 = pow(r0,1.0/3.0);
		r = r0/8.0;
79
		
80
		for(int k = 0; k < 64; ++k)
81
		{
82

83
84
85
			cnt[0] = cnt0[0] + 0.25*r0*sin(k/16.0*M_PI);
			cnt[1] = cnt0[1] + 0.25*r0*cos(k/16.0*M_PI);

86
			m.ClearFile();
87
88
89
90
91
92
			
			int numref;
			int refcnt = 0;
			do
			{
				numref = 0;
93
				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) 
94
				{
95
					if( am.GetLevel(it->self()) < max_levels )
96
					{
97
98
99
						it->Barycenter(xyz);
						//refine a circle
						q = 0;
100
						for(int d = 0; d < 2; ++d)
101
102
							q += (xyz[d]-cnt[d])*(xyz[d]-cnt[d]);
						q = sqrt(q);
103
104
						//if( q < r*(k+1) && q > r*k)
						if( q < r*1.4 && q > r)
105
106
107
108
						{
							indicator[*it] = 1;
							numref++;
						}
109
					}
110
					else indicator[*it] = 0;
111
				}
112
113
				numref = m.Integrate(numref);
				if( numref )
114
				{
115
116
					int ncells = m.TotalNumberOf(CELL);
					if( m.GetProcessorRank() == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
117
						std::cout << "k " << k << " refcnt " << refcnt << " " <<  r*k << " < r < " << r*(k+1) << " cells " << ncells << std::endl;
118
119
120
121
					//m.BeginSequentialCode();
					//std::cout << m.GetProcessorRank() << " cells " << m.NumberOfCells() << std::endl;
					//m.EndSequentialCode();

Kirill Terekhov's avatar
Kirill Terekhov committed
122
					if (!am.Refine(indicator)) break;
Kirill Terekhov's avatar
Kirill Terekhov committed
123
124
125
					
					if( false )
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
126
						std::stringstream file;
127
						file << "ref_" << k << "_" << refcnt << file_format;
Kirill Terekhov's avatar
Kirill Terekhov committed
128
129
						m.Save(file.str());
						if( m.GetProcessorRank() == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
130
							std::cout << "Save " << file.str() << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
131
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
132
					
133
134
135
136
137
138
139
140
				}
				refcnt++;
			}
			while(numref);
			refcnt = 0;
			do
			{
				numref = 0;
141
				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
142
				{
143
					if( am.GetLevel(it->self()) > 0 )
144
					{
145
146
147
						it->Barycenter(xyz);
						//refine a circle
						q = 0;
148
						for(int d = 0; d < 2; ++d)
149
150
							q += (xyz[d]-cnt[d])*(xyz[d]-cnt[d]);
						q = sqrt(q);
151
152
						//if( q < r*k)
						if( !(q < r*1.4 && q > r) )
153
154
155
156
						{
							indicator[*it] = 1;
							numref++;
						}
157
					}
158
					else indicator[*it] = 0;
159
				}
160
161
				numref = m.Integrate(numref);
				if( numref )
162
				{
163
164
					int ncells = m.TotalNumberOf(CELL);
					if( m.GetProcessorRank() == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
165
						std::cout << ": k " << k << " crscnt " << refcnt << " " << r*k << " < r < " << r*(k+1) << " cells " << ncells <<  std::endl;
166
167
168
					//m.BeginSequentialCode();
					//std::cout << m.GetProcessorRank() << " cells " << m.NumberOfCells() << std::endl;
					//m.EndSequentialCode();
Kirill Terekhov's avatar
Kirill Terekhov committed
169
170
					
					if( !am.Coarse(indicator) ) break;
Kirill Terekhov's avatar
Kirill Terekhov committed
171
172
173
174
					
					if( false ) 
					{
						std::stringstream file;
175
						file << "crs_" << k << "_" << refcnt << file_format;
Kirill Terekhov's avatar
Kirill Terekhov committed
176
177
178
179
180
						m.Save(file.str());
						if( m.GetProcessorRank() == 0 )
							std::cout << "Save " << file.str() << std::endl;
					}
					
181
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
182
				
183
184
185
186
				refcnt++;
			}
			while(numref);
			
Kirill Terekhov's avatar
Kirill Terekhov committed
187
			
188
#if defined(USE_PARTITIONER)
189
			if( true )
190
			{
191
192
				m.Barrier();
				std::cout << "before on " << m.GetProcessorRank() << " " << m.NumberOfCells() << std::endl;
193
				p.Evaluate();
194
				am.CheckParentSet(__FILE__,__LINE__);
195
196
				//am.Repartition();
				m.Redistribute();
197
198
199
				//std::fstream fout("sets"+std::to_string(m.GetProcessorRank())+".txt",std::ios::out);
				//am.ReportSets(fout);
				am.CheckParentSet(__FILE__,__LINE__);
200
				std::cout << "after on " << m.GetProcessorRank() << " " << m.NumberOfCells() << std::endl;
201
202
			}
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
203
			 
204
205
206
			 

			
Kirill Terekhov's avatar
Kirill Terekhov committed
207
			if( true )
208
			{
209
210
211
212
213
214
215
216
				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();
				}
217
				std::stringstream file;
218
				file << "step_" << k << file_format;
219
				m.Save(file.str());
Kirill Terekhov's avatar
Kirill Terekhov committed
220
221
				if( m.GetProcessorRank() == 0 )
					std::cout << "Save " << file.str() << std::endl;
222
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
223
			else if( m.GetProcessorRank() == 0 ) std::cout << "step " << k << std::endl;
224
225
226
		}
	}
	else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2]" << std::endl;
227
228
	
	Mesh::Finalize();
Kirill Terekhov's avatar
Kirill Terekhov committed
229
	return 0;
230
}