main.cpp 4.51 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
#include "amesh.h"

using namespace INMOST;

int main(int argc, char ** argv)
{
	Mesh::Initialize(&argc,&argv);
	
	if( argc > 1 )
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
11
		Mesh m;
12
		m.SetCommunicator(INMOST_MPI_COMM_WORLD);
Kirill Terekhov's avatar
Kirill Terekhov committed
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
		if( m.isParallelFileFormat(argv[1]) )
			m.Load(argv[1]);
		else if( m.GetProcessorRank() == 0 )
			m.Load(argv[1]);
	
#if defined(USE_PARTITIONER)
		if( true )
		{
			std::cout << "before on " << m.GetProcessorRank() << " " << m.NumberOfCells() << std::endl;
			Partitioner p(&m);
			p.SetMethod(Partitioner::INNER_KMEANS,Partitioner::Partition);
			p.Evaluate();
			m.Redistribute();
			m.ReorderEmpty(CELL|FACE|EDGE|NODE);
			std::cout << "after on " << m.GetProcessorRank() << " " << m.NumberOfCells() << std::endl;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
29
30
#endif
		m.ExchangeGhost(2,NODE);
Kirill Terekhov's avatar
Kirill Terekhov committed
31
		AdaptiveMesh am(m);
32
33
34
35
		//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
36
37
38
39
40
41
42
43
		/*
		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;
		*/
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
		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};
		Storage::real xyz[3], r, q, cnt[3];
		//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];
			}
		}
59
60
		m.AggregateMax(cmax,3);
		m.AggregateMin(cmin,3);
61
62
63
64
65
66
67
68
		r = 1;
		for(int d = 0; d < 3; ++d)
		{
			r *= cmax[d]-cmin[d];
			cnt[d] = (cmax[d]+cmin[d])*0.5;
		}
		r = pow(r,1.0/3.0)/20.0;
		
Kirill Terekhov's avatar
Kirill Terekhov committed
69
		for(int k = 0; k < 15; ++k)
70
71
72
73
74
75
76
		{
			
			int numref;
			int refcnt = 0;
			do
			{
				numref = 0;
77
				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) 
78
				{
79
					if( am.GetLevel(it->self()) < max_levels )
80
					{
81
82
83
84
85
86
87
88
89
90
91
						it->Barycenter(xyz);
						//refine a circle
						q = 0;
						for(int d = 0; d < 3; ++d)
							q += (xyz[d]-cnt[d])*(xyz[d]-cnt[d]);
						q = sqrt(q);
						if( q < r*(k+1) && q > r*k)
						{
							indicator[*it] = 1;
							numref++;
						}
92
					}
93
					else indicator[*it] = 0;
94
				}
95
				numref = m.Integrate(numref);
96
97
98
				if( numref )
				{
					std::cout << "k " << k << " refcnt " << refcnt << " " <<  r*k << " < r < " << r*(k+1) << " numref " << numref << " cells " << m.NumberOfCells() << std::endl;
99
100
					
					int res = am.Refine(indicator);
101
102
                    res = m.Integrate(res);
                    if (!res) break;
Kirill Terekhov's avatar
Kirill Terekhov committed
103
                    /*
Kirill Terekhov's avatar
Kirill Terekhov committed
104
105
106
107
108
109
110
                    {
						std::stringstream file;
						file << "ref_" << k << "_" << refcnt << ".pvtk";
						m.Save(file.str());
						if( m.GetProcessorRank() == 0 )
						std::cout << "Save " << file.str() << std::endl;
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
111
					 */
112
113
114
115
116
117
118
119
				}
				refcnt++;
			}
			while(numref);
			refcnt = 0;
			do
			{
				numref = 0;
120
				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
121
				{
122
					if( am.GetLevel(it->self()) > 0 )
123
					{
124
125
126
127
128
129
130
131
132
133
134
						it->Barycenter(xyz);
						//refine a circle
						q = 0;
						for(int d = 0; d < 3; ++d)
							q += (xyz[d]-cnt[d])*(xyz[d]-cnt[d]);
						q = sqrt(q);
						if( q < r*k)
						{
							indicator[*it] = 1;
							numref++;
						}
135
					}
136
					else indicator[*it] = 0;
137
				}
138
				numref = m.Integrate(numref);
139
140
				if( numref )
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
141
					std::cout << ": k " << k << " crscnt " << refcnt << " " << r*k << " < r < " << r*(k+1) << " numcrs " << numref << " cells " << m.NumberOfCells() <<  std::endl;
142
143
144
					int res = am.Coarse(indicator);
					res = m.Integrate(res);
					if( !res ) break;
145
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
146
				
147
148
149
150
151
152
				refcnt++;
			}
			while(numref);
			
			
			{
153
154
155
156
157
158
159
160
				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();
				}
161
				std::stringstream file;
Kirill Terekhov's avatar
Kirill Terekhov committed
162
				file << "step_" << k << ".pvtk";
163
				m.Save(file.str());
Kirill Terekhov's avatar
Kirill Terekhov committed
164
165
				if( m.GetProcessorRank() == 0 )
					std::cout << "Save " << file.str() << std::endl;
166
167
168
169
			}
		}
	}
	else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2]" << std::endl;
170
171
	
	Mesh::Finalize();
Kirill Terekhov's avatar
Kirill Terekhov committed
172
	return 0;
173
}