amesh.cpp 46.9 KB
Newer Older
1
#include "amesh.h"
SilverLife's avatar
SilverLife committed
2
#include <iomanip>
3
#include <set>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25


#if defined(USE_PARALLEL_WRITE_TIME)
#define REPORT_MPI(x) {m->WriteTab(m->out_time) << "<MPI><![CDATA[" << #x << "]]></MPI>" << std::endl; x;}
#define REPORT_STR(x) {m->WriteTab(m->out_time) << "<TEXT><![CDATA[" << x << "]]></TEXT>" << std::endl;}
#define REPORT_VAL(str,x) {m->WriteTab(m->out_time) << "<VALUE name=\"" << str << "\"> <CONTENT><![CDATA[" << x << "]]></CONTENT> <CODE><![CDATA[" << #x << "]]></CODE></VALUE>" << std::endl;}
#define ENTER_FUNC() double all_time = Timer(); {m->WriteTab(m->out_time) << "<FUNCTION name=\"" << __FUNCTION__ << "\" id=\"func" << m->GetFuncID()++ << "\">" << std::endl; m->Enter();}
#define ENTER_BLOCK() { double btime = Timer(); m->WriteTab(m->out_time) << "<FUNCTION name=\"" << __FUNCTION__ << ":" << __FILE__ << ":" << __LINE__ << "\" id=\"func" << m->GetFuncID()++ << "\">" << std::endl; m->Enter();
#define EXIT_BLOCK() m->WriteTab(m->out_time) << "<TIME>" << Timer() - btime << "</TIME>" << std::endl; m->Exit(); m->WriteTab(m->out_time) << "</FUNCTION>" << std::endl;}
#define EXIT_FUNC() {m->WriteTab(m->out_time) << "<TIME>" << Timer() - all_time << "</TIME>" << std::endl; m->Exit(); m->WriteTab(m->out_time) << "</FUNCTION>" << std::endl;}
#define EXIT_FUNC_DIE() {m->WriteTab(m->out_time) << "<TIME>" << -1 << "</TIME>" << std::endl; m->Exit(); m->WriteTab(m->out_time) << "</FUNCTION>" << std::endl;}
#else
#define REPORT_MPI(x) x
#define REPORT_STR(x) {}
#define REPORT_VAL(str,x) {}
#define ENTER_FUNC() {}
#define EXIT_FUNC() {}
#define EXIT_FUNC_DIE()  {}
#endif


//#include "../../Source/Misc/base64.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
26
//using namespace std;
27

Kirill Terekhov's avatar
Kirill Terekhov committed
28
29
//from inmost
//std::string base64_encode(unsigned char const* buf, unsigned int bufLen);
Kirill Terekhov's avatar
Kirill Terekhov committed
30
/*
Kirill Terekhov's avatar
Kirill Terekhov committed
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
	std::string base64_encode_(unsigned char const* buf, unsigned int bufLen)
	{
		static const std::string base64_chars =
		"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
		"abcdefghijklmnopqrstuvwxyz"
		"0123456789+/";
		std::string ret;
		int i = 0;
		int j = 0;
		unsigned char char_array_3[3];
		unsigned char char_array_4[4];
		
		while (bufLen--)
		{
			char_array_3[i++] = *(buf++);
			if (i == 3)
			{
				char_array_4[0] = (char_array_3[0] & 0xfc) >> 2;
				char_array_4[1] = ((char_array_3[0] & 0x03) << 4) + ((char_array_3[1] & 0xf0) >> 4);
				char_array_4[2] = ((char_array_3[1] & 0x0f) << 2) + ((char_array_3[2] & 0xc0) >> 6);
				char_array_4[3] = char_array_3[2] & 0x3f;
				for(i = 0; (i <4) ; i++)
					ret += base64_chars[char_array_4[i]];
				i = 0;
			}
		}
		
		if (i)
		{
			for(j = i; j < 3; j++)
				char_array_3[j] = '\0';
	  
			char_array_4[0] = (char_array_3[0] & 0xfc) >> 2;
			char_array_4[1] = ((char_array_3[0] & 0x03) << 4) + ((char_array_3[1] & 0xf0) >> 4);
			char_array_4[2] = ((char_array_3[1] & 0x0f) << 2) + ((char_array_3[2] & 0xc0) >> 6);
			char_array_4[3] = char_array_3[2] & 0x3f;
	  
			for (j = 0; (j < i + 1); j++)
				ret += base64_chars[char_array_4[j]];
	  
			while((i++ < 3))
				ret += '=';
		}
		
		return ret;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
77
 */
78
79
80
/// todo:
/// 1. coarsment
/// 2. strategy for faces/edges with faults
Kirill Terekhov's avatar
Kirill Terekhov committed
81
/// 3. geom model supportbn
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
/// 4. make class abstract virtual for user implementation of refinement and coarsment indicators
/// see in code todo:
namespace INMOST
{
	void CleanupSets(ElementSet set)
	{
		ElementSet::iterator it = set.Begin();
		while(it != set.End())
		{
			if( it->isValid() ) ++it;
			else it = set.Erase(it);
		}
		for(ElementSet child = set.GetChild(); child.isValid(); child = child.GetSibling())
			CleanupSets(child);
	}
	
	void ReduceMax(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void) size;
		element->Integer(tag) = std::max(element->Integer(tag),*((const INMOST_DATA_INTEGER_TYPE *)data));
	}
	
	void ReduceMin(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void) size;
		element->Integer(tag) = std::min(element->Integer(tag),*((const INMOST_DATA_INTEGER_TYPE *)data));
	}
SilverLife's avatar
SilverLife committed
109

110
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
111
    void AdaptiveMesh::PrintSetLocal(std::string offset, ElementSet it, std::stringstream& ss)
SilverLife's avatar
SilverLife committed
112
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
113
114
        std::stringstream ss1;
        ss1 << offset << rank << ": Set : " << std::setw(5) << it.GetName() << " ";
SilverLife's avatar
SilverLife committed
115
116
        for (int i = ss1.str().length(); i < 23; i++) ss1 << " ";
        ss << ss1.str();
Kirill Terekhov's avatar
Kirill Terekhov committed
117
        ss << std::setw(6);
SilverLife's avatar
SilverLife committed
118
119
120
121
122
        if (it.GetStatus() == Element::Shared) ss << "shared";
        else if (it.GetStatus() == Element::Ghost) ss << "ghost";
        else if (it.GetStatus() == Element::Owned) ss << "owned";
        else ss << "none";

123
        ss << " tag_owner (" << it.IntegerDF(m->OwnerTag()) << ")";
SilverLife's avatar
SilverLife committed
124
125
126

        //ss << "   level (" << level[it.self()] << ")  ";
        ss << " tag_processors (";
Kirill Terekhov's avatar
Kirill Terekhov committed
127
        std::stringstream ss2;
128
        Storage::integer_array arr = it.IntegerArrayDV(m->ProcessorsTag());
SilverLife's avatar
SilverLife committed
129
130
        for (int i = 0; i < arr.size(); i++)
            ss2 << arr[i] << " ";
Kirill Terekhov's avatar
Kirill Terekhov committed
131
        ss << std::setw(5) << ss2.str() <<")";
SilverLife's avatar
SilverLife committed
132
133
134
135
136
137
138
139
        

        ElementSet::iterator p = it.Begin();
        ss << "     | Refs: ";
        int first = 0;
        while(p != it.End())
        {
        //    if (first++ == 0) ss << endl << offset << "   ";
Kirill Terekhov's avatar
Kirill Terekhov committed
140
            std::string type = "unknw";
SilverLife's avatar
SilverLife committed
141
142
143
144
            if (p->GetElementType() == CELL) type = "cell";
            if (p->GetElementType() == FACE) type = "face";
            if (p->GetElementType() == EDGE) type = "edge";
            if (p->GetElementType() == NODE) type = "node";
Kirill Terekhov's avatar
Kirill Terekhov committed
145
            ss << type << "-" << std::setw(2) << p->GlobalID() << " ";
SilverLife's avatar
SilverLife committed
146
147
            p++;
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
148
        ss << std::endl;
SilverLife's avatar
SilverLife committed
149
150
151
152
153
154

        for(ElementSet child = it.GetChild(); child.isValid(); child = child.GetSibling())
        {
            PrintSetLocal(offset + "   ",child,ss);
        }
    }
155
	*/
SilverLife's avatar
SilverLife committed
156

157
	/*
SilverLife's avatar
SilverLife committed
158
159
    void AdaptiveMesh::PrintSet()
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
160
        std::stringstream ss;
161
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
SilverLife's avatar
SilverLife committed
162
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
163
            //if (it->HaveParent()) continue;
164
            PrintSetLocal("",ElementSet(m,*it),ss);
SilverLife's avatar
SilverLife committed
165
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
166
        std::cout << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
167
    }
168
	*/
SilverLife's avatar
SilverLife committed
169

170
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
171
    void PrintRefs(std::ostream& os, Storage::reference_array refs)
172
173
174
    {
        for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
175
            std::string type = "unknw";
176
177
178
179
180
181
182
            if (refs[i].GetElementType() == CELL) type = "cell";
            if (refs[i].GetElementType() == FACE) type = "face";
            if (refs[i].GetElementType() == EDGE) type = "edge";
            if (refs[i].GetElementType() == NODE) type = "node";
            os << "(" << type << "," << refs[i]->GlobalID() << ") ";
			Storage::real xyz[3] = {0,0,0};
            refs[i]->Centroid(xyz);
Kirill Terekhov's avatar
Kirill Terekhov committed
183
            os << "(" << xyz[0] << "," << xyz[1] << "," << xyz[2] <<")" <<  std::endl;
184
185
        }
    }
186
	*/
187

188
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
189
    void AdaptiveMesh::PrintMesh(std::ostream& os, int cell, int face, int edge, int node)
SilverLife's avatar
SilverLife committed
190
191
    {
        if (cell + face + edge + node == 0) return;
Kirill Terekhov's avatar
Kirill Terekhov committed
192
193
        std::stringstream ss;
        ss << "================= " << rank << " =====================" << std::endl;
SilverLife's avatar
SilverLife committed
194
195
        if (cell)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
196
            ss << "Cells: " << m->NumberOfCells() <<  std::endl;
197
            for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) 
SilverLife's avatar
SilverLife committed
198
            {
199
                ss << rank << ": " << it->GlobalID() << " - " << it->LocalID() << " - ";            
SilverLife's avatar
SilverLife committed
200
201
202
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";
Kirill Terekhov's avatar
Kirill Terekhov committed
203
                ss << std::endl;
SilverLife's avatar
SilverLife committed
204
205
206
207
208
            }
        }

        if (face)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
209
            ss << "Faces: " << m->NumberOfFaces() <<  std::endl;
210
            for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) 
SilverLife's avatar
SilverLife committed
211
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
212
213
                ss << rank << ": " << std::setw(2) << it->LocalID() << " " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
214
215
216
217
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";

218
219
220

                double xyz[3];
                it->Centroid(xyz);
Kirill Terekhov's avatar
Kirill Terekhov committed
221
                ss << "   (" << std::setw(5) << xyz[0] << " " << std::setw(5) << xyz[1] << " " << std::setw(5) << xyz[2] << ")";
222

223
                ss << "  " << m->GetMarker(*it,m->NewMarker());
224
225
226
227

                ss << " nc(" << it->getNodes().size() << ": "; 
                ElementArray<Node> nodes = it->getNodes();
                for (ElementArray<Node>::iterator node = nodes.begin(); node != nodes.end(); node++)
Kirill Terekhov's avatar
Kirill Terekhov committed
228
                    ss << std::setw(2) << node->GlobalID() << " ";
229
                ss << ")";
Kirill Terekhov's avatar
Kirill Terekhov committed
230
                ss << std::endl;
SilverLife's avatar
SilverLife committed
231
232
233
234
235
            }
        }

        if (edge)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
236
            ss << "Edges: " << m->NumberOfEdges() <<  std::endl;
237
            for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it) 
SilverLife's avatar
SilverLife committed
238
239
240
241
242
            {
                ss << rank << ": " << it->GlobalID() << " - " ;            
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";
Kirill Terekhov's avatar
Kirill Terekhov committed
243
                ss << std::endl;
SilverLife's avatar
SilverLife committed
244
245
246
247
248
            }
        }

        if (node)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
249
            ss << "Nodes:" << std::endl;
250
            for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) 
SilverLife's avatar
SilverLife committed
251
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
252
253
                ss << rank << ": " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
254
255
256
257
258
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";

                {
259
                    ss << "  " << m->GetMarker(*it,m->NewMarker());
260
261

                    ss << "(" << 
Kirill Terekhov's avatar
Kirill Terekhov committed
262
263
264
                        std::setw(3) << it->RealArray(m->CoordsTag())[0] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[1] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[2] << ")";
265

SilverLife's avatar
SilverLife committed
266
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
267
                ss << std::endl;
SilverLife's avatar
SilverLife committed
268
269
270
            }
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
271
272
        ss << "=========================================" << std::endl;
        os << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
273
    }
274
	*/
SilverLife's avatar
SilverLife committed
275

276
	/*
SilverLife's avatar
SilverLife committed
277
278
279
280
281
    void AdaptiveMesh::UpdateStatus()
    {

        for(ElementType mask = CELL; mask >= NODE; mask = PrevElementType(mask))
        {
282
            for(Mesh::iteratorElement it = m->BeginElement(mask); it != m->EndElement(); it++)
SilverLife's avatar
SilverLife committed
283
284
285
286
287
288
289
290
291
            {
                int stat = 0;
                if (it->GetStatus() == Element::Shared) stat = 1;
                else if (it->GetStatus() == Element::Ghost)  stat = 2;

                tag_status[it->self()] = stat;
            }
        }
    }
292
	*/
SilverLife's avatar
SilverLife committed
293
    
294
	/*
SilverLife's avatar
SilverLife committed
295
296
297
298
299
300
301
    void AdaptiveMesh::PrintSet(ElementSet set, std::string offset)
    {
        std::cout << offset << "Set: " << std::endl;

		ElementSet::iterator it = set.Begin();
		while(it != set.End())
        {
302
		    ElementSet parent(m,parent_set[*it]);
SilverLife's avatar
SilverLife committed
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
            std::cout << offset << it->GlobalID() << " - " << level[*it] << " : ";
		    
            ElementSet::iterator p = parent.Begin();
    		while(p != parent.End())
            {
                std::cout << p->GlobalID() << " ";
                p++;
            }
            std::cout << std::endl;
            it++;
        }

		for(ElementSet child = set.GetChild(); child.isValid(); child = child.GetSibling())
        {
            PrintSet(child,offset + "  ");
        }

    }
321
	*/
SilverLife's avatar
SilverLife committed
322

323
	/*
324
325
326
    void AdaptiveMesh::SynchronizeSet(ElementSet set)
    {
    #ifdef USE_MPI
327
328
        int size = m->GetProcessorsNumber();
        int rank = m->GetProcessorRank();
329
330
        for (int i = 0; i < size; i++)
        {
331
332
            set.IntegerArray(m->SendtoTag()).push_back(i);
            m->ExchangeMarked();
333
334
335
        }
    #endif
    }
336
	*/
337

338
	/*
SilverLife's avatar
SilverLife committed
339
340
    void AdaptiveMesh::Test()
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
341
        std::cout << rank << ": ================" << std::endl;
SilverLife's avatar
SilverLife committed
342
343
        PrintSet(root,"");
    }
344
	*/
345
346
347
	
	void AdaptiveMesh::ClearData()
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
348
		set_id        = m->DeleteTag(set_id);
Kirill Terekhov's avatar
Kirill Terekhov committed
349
350
351
		level         = m->DeleteTag(level);
		hanging_nodes = m->DeleteTag(hanging_nodes);
		parent_set    = m->DeleteTag(parent_set);
352
353
354
355
356
357
358
359
		root.DeleteSetTree();
	}
	
	void AdaptiveMesh::PrepareSet()
	{
		//retrive set for coarsening, initialize set if is not present
		if( !root.isValid() )
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
360
			root = m->GetSet("ROOT_SET");
361
362
			if( root == InvalidElement() )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
363
				root = m->CreateSetUnique("ROOT_SET").first;
Kirill Terekhov's avatar
Kirill Terekhov committed
364
				root.SetExchange(ElementSet::SYNC_ELEMENTS_SHARED);
365
				level[root] = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
366
				for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
367
368
369
370
				{
					root.PutElement(it->self());
					parent_set[it->self()] = root.GetHandle();
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
371
				m->Enumerate(CELL,set_id);
372
			}
SilverLife's avatar
SilverLife committed
373
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
374
		if( !m->HaveGlobalID(CELL) ) m->AssignGlobalID(CELL); //for unique set names
Kirill Terekhov's avatar
Kirill Terekhov committed
375
		m->ResolveSets();
376
377
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
378
	AdaptiveMesh::AdaptiveMesh(Mesh & _m) : m(&_m)
379
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
380
		model = NULL;
381
		//create a tag that stores maximal refinement level of each element
Kirill Terekhov's avatar
Kirill Terekhov committed
382
		level = m->CreateTag("REFINEMENT_LEVEL",DATA_INTEGER,CELL|FACE|EDGE|NODE|ESET,NONE,1);
383
		//tag_status = m->CreateTag("TAG_STATUS",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
384
		set_id = m->CreateTag("SET_ID",DATA_INTEGER,CELL,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
385
386
		//tag_an = m->CreateTag("TAG_AN",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
		//ref_tag = m->CreateTag("REF",DATA_REFERENCE,CELL|FACE|EDGE|NODE,NONE);
387
		//create a tag that stores links to all the hanging nodes of the cell
Kirill Terekhov's avatar
Kirill Terekhov committed
388
		hanging_nodes = m->CreateTag("HANGING_NODES",DATA_REFERENCE,CELL|FACE,NONE);
389
		//create a tag that stores links to sets
Kirill Terekhov's avatar
Kirill Terekhov committed
390
		parent_set = m->CreateTag("PARENT_SET",DATA_REFERENCE,CELL,NONE,1);
391
392
	    size = m->GetProcessorsNumber();
    	rank = m->GetProcessorRank();
393
394
395
396
397
398
399
400
	}
	
	AdaptiveMesh::~AdaptiveMesh()
	{
		//do not delete tags, user may want to repetitively use this class
		//as extension of class mesh in limited code span
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
401
402
	void AdaptiveMesh::CheckParentSet()
	{
403
		ENTER_FUNC();
Kirill Terekhov's avatar
Kirill Terekhov committed
404
405
406
407
408
		int err = 0;
		for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
		{
			if( parent_set[*it] == InvalidHandle() )
			{
409
				REPORT_STR(m->GetProcessorRank() << " parent set not valid on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << " " << level[*it]);
Kirill Terekhov's avatar
Kirill Terekhov committed
410
411
412
413
				err++;
			}
			else if( GetHandleElementType(parent_set[*it]) != ESET )
			{
414
				REPORT_STR(m->GetProcessorRank() << " parent set is something else " << ElementTypeName(GetHandleElementType(parent_set[*it])) << ":" << GetHandleID(parent_set[*it]) << " on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << " " << level[*it]);
Kirill Terekhov's avatar
Kirill Terekhov committed
415
416
417
418
				err++;
			}
		}
		err = m->Integrate(err);
419
		EXIT_FUNC();
Kirill Terekhov's avatar
Kirill Terekhov committed
420
421
		if( err ) 
		{
422
			REPORT_STR(rank << " error in " << __FUNCTION__);
Kirill Terekhov's avatar
Kirill Terekhov committed
423
424
425
426
427
428
			std::cout << rank << " error in " << __FUNCTION__ << std::endl;
			exit(-1);
		}
		
	}
	
429
430
	bool AdaptiveMesh::Refine(TagInteger & indicator)
	{
431
432
		static int fi = 0;
        ENTER_FUNC();
433
434
435
		static int call_counter = 0;
		int ret = 0; //return number of refined cells
		//initialize tree structure
436
437
		//m->CheckCentroids(__FILE__,__LINE__);
		ENTER_BLOCK();
438
		PrepareSet();
439
440
441
442
		EXIT_BLOCK();

		m->Save("before_refine"+std::to_string(fi)+".pvtk");
		std::cout << "Save before_refine"+std::to_string(fi)+".pvtk" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
443
		
444
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
445
446
447
448
449
450
451
452
453
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(parent_set,CELL,0);
		m->CheckCentroids(__FILE__,__LINE__);
		m->ResolveSets();
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(parent_set,CELL,0);
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
		m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
454
		CheckParentSet();
455
456
457
458
459
460
461
		EXIT_BLOCK();

		m->Save("before_refine_parent"+std::to_string(fi)+".pvtk");
		std::cout << "Save before_refine_parent"+std::to_string(fi)+".pvtk" << std::endl;
		

		
462
463
		int schedule_counter = 1; //indicates order in which refinement will be scheduled
		int scheduled = 1; //indicates that at least one element was scheduled on current sweep
464
		ENTER_BLOCK();
465
		//0. Extend indicator for edges and faces
Kirill Terekhov's avatar
Kirill Terekhov committed
466
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
467
468
469
		while(scheduled)
		{
			//1.Communicate indicator - it may be not synced
Kirill Terekhov's avatar
Kirill Terekhov committed
470
			m->ExchangeData(indicator,CELL,0);
471
472
			//2.Propogate indicator down to the faces,edges
			//  select schedule for them
Kirill Terekhov's avatar
Kirill Terekhov committed
473
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
474
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
475
				Cell c = m->CellByLocalID(it);
476
477
478
479
480
481
482
483
484
485
486
				if( indicator[c] == schedule_counter )
				{
					ElementArray<Element> adj = c.getAdjElements(FACE|EDGE);
					for(ElementArray<Element>::size_type kt = 0; kt < adj.size(); ++kt)
					{
						if( level[adj[kt]] == level[c] ) //do not schedule finer or coarser elements
							indicator[adj[kt]] = schedule_counter; //refine together with current cell
					}
				}
			}
			//3.Communicate indicator on faces and edges
Kirill Terekhov's avatar
Kirill Terekhov committed
487
			m->ExchangeData(indicator,FACE|EDGE,0);
488
489
490
491
			//4.Check for each cell if there is
			//  any hanging node with adjacent in a need to refine,
			//  schedule for refinement earlier.
			scheduled = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
492
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
493
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
494
				Cell c = m->CellByLocalID(it);
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
				//already scheduled cells may be required to be refined first
				//if( indicator[c] == 0 ) //some optimization
				{
					bool scheduled_c = false;
					//any finer level edge is scheduled to be refined first
					ElementArray<Edge> edges = c->getEdges();
					for(ElementArray<Edge>::size_type kt = 0; kt < edges.size() && !scheduled_c; ++kt)
					{
						//if a finer edge is scheduled
						//then this cell should be refined first
						if( indicator[edges[kt]] != 0 &&
							level[edges[kt]] > level[c] &&
							indicator[edges[kt]] >= indicator[c] )
						{
							indicator[c] = schedule_counter+1;
							scheduled++;
							scheduled_c = true;
						}
					}
				}
			}
			//5.Go back to 1 until no new elements scheduled
Kirill Terekhov's avatar
Kirill Terekhov committed
517
			scheduled = m->Integrate(scheduled);
518
519
			if( scheduled ) schedule_counter++;
		}
520
		m->ExchangeData(indicator,CELL | FACE | EDGE,0);
521
		EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
522
		//m->Save("indicator.pmf");
523
		//6.Refine
524
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
525
		m->BeginModification();
526
527
528
529
		while(schedule_counter)
		{
			Storage::real xyz[3] = {0,0,0};
			//7.split all edges of the current schedule
Kirill Terekhov's avatar
Kirill Terekhov committed
530
			for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
531
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
532
				Edge e = m->EdgeByLocalID(it);
533
534
535
536
537
				if( !e.Hidden() && indicator[e] == schedule_counter )
				{
					//remember adjacent faces that should get information about new hanging node
					ElementArray<Face> edge_faces = e.getFaces();
					//location on the center of the edge
Kirill Terekhov's avatar
Kirill Terekhov committed
538
					for(Storage::integer d = 0; d < m->GetDimensions(); ++d)
539
540
541
						xyz[d] = (e.getBeg().Coords()[d]+e.getEnd().Coords()[d])*0.5;
					//todo: request transformation of node location according to geometrical model
					//create middle node
Kirill Terekhov's avatar
Kirill Terekhov committed
542
					Node n = m->CreateNode(xyz);
543
					//set increased level for new node
544
					level[n] = level[e.getBeg()] = level[e.getEnd()] = level[e]+1;
545
546
547
548
					//for each face provide link to a new hanging node
					for(ElementArray<Face>::size_type kt = 0; kt < edge_faces.size(); ++kt)
						hanging_nodes[edge_faces[kt]].push_back(n);
					//split the edge by the middle node
Kirill Terekhov's avatar
Kirill Terekhov committed
549
					ElementArray<Edge> new_edges = Edge::SplitEdge(e,ElementArray<Node>(m,1,n.GetHandle()),0);
550
551
552
553
554
					//set increased level for new edges
					level[new_edges[0]] = level[new_edges[1]] = level[e]+1;
				}
			}
			//8.split all faces of the current schedule, using hanging nodes on edges
Kirill Terekhov's avatar
Kirill Terekhov committed
555
			for(Storage::integer it = 0; it < m->FaceLastLocalID(); ++it) if( m->isValidFace(it) )
556
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
557
				Face f = m->FaceByLocalID(it);
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
				if( !f.Hidden() && indicator[f] == schedule_counter )
				{
					//connect face center to hanging nodes of the face
					Storage::reference_array face_hanging_nodes = hanging_nodes[f];
					//remember adjacent cells that should get information about new hanging node
					//and new hanging edges
					ElementArray<Cell> face_cells = f.getCells();
					//create node at face center
					//f->Centroid(xyz);
					for(int d = 0; d < 3; ++d) xyz[d] = 0.0;
					for(Storage::reference_array::size_type kt = 0; kt < face_hanging_nodes.size(); ++kt)
						for(int d = 0; d < 3; ++d) xyz[d] += face_hanging_nodes[kt].getAsNode().Coords()[d];
					for(int d = 0; d < 3; ++d) xyz[d] /= (Storage::real)face_hanging_nodes.size();
					//todo: request transformation of node location according to geometrical model
					//create middle node
Kirill Terekhov's avatar
Kirill Terekhov committed
573
					Node n = m->CreateNode(xyz);
574
					//set increased level for the new node
575
					level[n] = level[f]+1;
576
577
578
					//for each cell provide link to new hanging node
					for(ElementArray<Face>::size_type kt = 0; kt < face_cells.size(); ++kt)
						hanging_nodes[face_cells[kt]].push_back(n);
Kirill Terekhov's avatar
Kirill Terekhov committed
579
580
					ElementArray<Node> edge_nodes(m,2); //to create new edges
					ElementArray<Edge> hanging_edges(m,face_hanging_nodes.size());
581
582
583
584
					edge_nodes[0] = n;
					for(Storage::reference_array::size_type kt = 0; kt < face_hanging_nodes.size(); ++kt)
					{
						edge_nodes[1] = face_hanging_nodes[kt].getAsNode();
Kirill Terekhov's avatar
Kirill Terekhov committed
585
						hanging_edges[kt] = m->CreateEdge(edge_nodes).first;
586
587
588
589
590
591
592
593
594
595
596
						//set increased level for new edges
						level[hanging_edges[kt]] = level[f]+1;
					}
					//split the face by these edges
					ElementArray<Face> new_faces = Face::SplitFace(f,hanging_edges,0);
					//set increased level to new faces
					for(ElementArray<Face>::size_type kt = 0; kt < new_faces.size(); ++kt)
						level[new_faces[kt]] = level[f]+1;
				}
			}
			//this tag helps recreate internal face
Kirill Terekhov's avatar
Kirill Terekhov committed
597
			TagReferenceArray internal_face_edges = m->CreateTag("INTERNAL_FACE_EDGES",DATA_REFERENCE,NODE,NODE,4);
598
			//this marker helps detect edges of current cell only
Kirill Terekhov's avatar
Kirill Terekhov committed
599
			MarkerType mark_cell_edges = m->CreateMarker();
600
			//this marker helps detect nodes hanging on edges of unrefined cell
Kirill Terekhov's avatar
Kirill Terekhov committed
601
			MarkerType mark_hanging_nodes = m->CreateMarker();
602
			//9.split all cells of the current schedule
Kirill Terekhov's avatar
Kirill Terekhov committed
603
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
604
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
605
				Cell c = m->CellByLocalID(it);
606
607
608
609
610
611
612
613
614
615
616
				if( !c.Hidden() && indicator[c] == schedule_counter )
				{
					Storage::reference_array cell_hanging_nodes = hanging_nodes[c]; //nodes to be connected
					//create node at cell center
					for(int d = 0; d < 3; ++d) xyz[d] = 0.0;
					for(Storage::reference_array::size_type kt = 0; kt < cell_hanging_nodes.size(); ++kt)
						for(int d = 0; d < 3; ++d) xyz[d] += cell_hanging_nodes[kt].getAsNode().Coords()[d];
					for(int d = 0; d < 3; ++d) xyz[d] /= (Storage::real)cell_hanging_nodes.size();
					//c->Centroid(xyz);
					//todo: request transformation of node location according to geometrical model
					//create middle node
Kirill Terekhov's avatar
Kirill Terekhov committed
617
					Node n = m->CreateNode(xyz);
618
					//set increased level for the new node
619
					level[n] = level[c]+1;
620
621
622
623
624
					//retrive all edges of current face to mark them
					ElementArray<Edge> cell_edges = c.getEdges();
					//mark all edges so that we can retive them later
					cell_edges.SetMarker(mark_cell_edges);
					//connect face center to centers of faces by edges
Kirill Terekhov's avatar
Kirill Terekhov committed
625
626
					ElementArray<Node> edge_nodes(m,2);
					ElementArray<Edge> edges_to_faces(m,cell_hanging_nodes.size());
627
628
629
630
631
632
					edge_nodes[0] = n;
					for(Storage::reference_array::size_type kt = 0; kt < cell_hanging_nodes.size(); ++kt)
					{
						assert(cell_hanging_nodes[kt].isValid());
						//todo: unmark hanging node on edge if no more cells depend on it
						edge_nodes[1] = cell_hanging_nodes[kt].getAsNode();
Kirill Terekhov's avatar
Kirill Terekhov committed
633
						edges_to_faces[kt] = m->CreateEdge(edge_nodes).first;
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
						//set increased level for new edges
						level[edges_to_faces[kt]] = level[c]+1;
						//for each node other then the hanging node of the face
						//(this is hanging node on the edge)
						//we record a pair of edges to reconstruct internal faces
						ElementArray<Edge> hanging_edges = cell_hanging_nodes[kt].getEdges(mark_cell_edges,0);
						for(ElementArray<Edge>::size_type lt = 0; lt < hanging_edges.size(); ++lt)
						{
							//get hanging node on the edge
							assert(hanging_edges[lt].getBeg() == cell_hanging_nodes[kt] || hanging_edges[lt].getEnd() == cell_hanging_nodes[kt]);
							Node v = hanging_edges[lt].getBeg() == cell_hanging_nodes[kt]? hanging_edges[lt].getEnd() : hanging_edges[lt].getBeg();
							//mark so that we can collect all of them
							v.SetMarker(mark_hanging_nodes);
							//fill the edges
							Storage::reference_array face_edges = internal_face_edges[v];
							//fill first two in forward order
							//this way we make a closed loop
							assert(face_edges[0] == InvalidElement() || face_edges[2] == InvalidElement());
							if( face_edges[0] == InvalidElement() )
							{
								face_edges[0] = edges_to_faces[kt];
								face_edges[1] = hanging_edges[lt];
							}
							else //last two in reverse
							{
								assert(face_edges[2] ==InvalidElement());
								face_edges[2] = hanging_edges[lt];
								face_edges[3] = edges_to_faces[kt];
							}
						}
					}
					//remove marker from cell edges
					cell_edges.RemMarker(mark_cell_edges);
					//now we have to create internal faces
					ElementArray<Node> edge_hanging_nodes = c.getNodes(mark_hanging_nodes,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
669
					ElementArray<Face> internal_faces(m,edge_hanging_nodes.size());
670
671
672
673
674
675
676
677
678
679
					//unmark hanging nodes on edges
					edge_hanging_nodes.RemMarker(mark_hanging_nodes);
					for(ElementArray<Node>::size_type kt = 0; kt < edge_hanging_nodes.size(); ++kt)
					{
						//create a face based on collected edges
						Storage::reference_array face_edges = internal_face_edges[edge_hanging_nodes[kt]];
						assert(face_edges[0].isValid());
						assert(face_edges[1].isValid());
						assert(face_edges[2].isValid());
						assert(face_edges[3].isValid());
Kirill Terekhov's avatar
Kirill Terekhov committed
680
						internal_faces[kt] = m->CreateFace(ElementArray<Edge>(m,face_edges.begin(),face_edges.end())).first;
681
682
683
684
685
						//set increased level
						level[internal_faces[kt]] = level[c]+1;
						//clean up structure, so that other cells can use it
						edge_hanging_nodes[kt].DelData(internal_face_edges);
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
686
687
688
689
690
691
692
					//if( c.GlobalID() == 228 )
					//{
					//	double cnt[3];
					//	c.Centroid(cnt);
					//	std::cout << "Split CELL:" << c.LocalID() << " " << c.GlobalID() << " " << Element::StatusName(c.GetStatus()) << " " << cnt[0] << " " << cnt[1] << " " << cnt[2] << std::endl;
					//	
					//}
693
694
695
					//split the cell
					ElementArray<Cell> new_cells = Cell::SplitCell(c,internal_faces,0);
					//retrive parent set
Kirill Terekhov's avatar
Kirill Terekhov committed
696
					ElementSet parent(m,parent_set[c]);
697
					//create set corresponding to old coarse cell
Kirill Terekhov's avatar
Kirill Terekhov committed
698
699
					Storage::real cnt[3];
					c.Centroid(cnt);
700
					std::stringstream set_name;
Kirill Terekhov's avatar
Kirill Terekhov committed
701
					//set_name << parent.GetName() << "_C" << c.GlobalID(); //rand may be unsafe
Kirill Terekhov's avatar
Kirill Terekhov committed
702
703
704
705
					if( parent == root )
						set_name << "r" << set_id[c];
					else
						set_name << parent.GetName() << "c" << set_id[c];
706
					//set_name << base64_encode_((unsigned char *)cnt,3*sizeof(double)/sizeof(unsigned char));
Kirill Terekhov's avatar
Kirill Terekhov committed
707
					
Kirill Terekhov's avatar
Kirill Terekhov committed
708
709
710
711
712
713
714
715
716
717
					ElementSet check_set = m->GetSet(set_name.str());
					if( check_set.isValid() )
					{
						std::cout << rank << " set " << set_name.str() << " for cell " << c.GlobalID() << " " << Element::StatusName(c.GetStatus()) << " already exists" << std::endl;
						if( check_set->HaveParent() )
							std::cout << rank << " parent is " << check_set->GetParent()->GetName() << " cell parent is " << parent.GetName() << std::endl;
						std::cout << rank << " Elements: ";
						for(ElementSet::iterator it = check_set.Begin(); it != check_set.End(); ++it)
							std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << " ";
						std::cout << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
718
						exit(-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
719
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
720
					
Kirill Terekhov's avatar
Kirill Terekhov committed
721
					ElementSet cell_set = m->CreateSetUnique(set_name.str()).first;
Kirill Terekhov's avatar
Kirill Terekhov committed
722
					cell_set->SetExchange(ElementSet::SYNC_ELEMENTS_ALL);
723
724
725
726
					level[cell_set] = level[c]+1;
					//set up increased level for the new cells
					for(ElementArray<Cell>::size_type kt = 0; kt < new_cells.size(); ++kt)
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
727
						set_id[new_cells[kt]] = kt;
728
729
730
731
						level[new_cells[kt]] = level[c]+1;
						cell_set.PutElement(new_cells[kt]);
						parent_set[new_cells[kt]] = cell_set.GetHandle();
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
732
733
734
735
736
737
738
739
740
741
					/*
					if( check_set.isValid() )
					{
						std::cout << rank << " Elements: ";
						for(ElementSet::iterator it = check_set.Begin(); it != check_set.End(); ++it)
							std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << " ";
						std::cout << std::endl;
					}
					*/
					//if( !cell_set->HaveParent() )
742
					parent.AddChild(cell_set);
Kirill Terekhov's avatar
Kirill Terekhov committed
743
					//else assert(cell_set->GetParent() == parent);
744
745
746
747
					//increment number of refined cells
					ret++;
				}
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
748
749
750
			m->ReleaseMarker(mark_hanging_nodes);
			m->ReleaseMarker(mark_cell_edges);
			m->DeleteTag(internal_face_edges);
751
752
753
			//10.jump to later schedule, and go to 7.
			schedule_counter--;
		}
SilverLife's avatar
SilverLife committed
754

755
		//free created tag
Kirill Terekhov's avatar
Kirill Terekhov committed
756
		m->DeleteTag(indicator,FACE|EDGE);
SilverLife's avatar
SilverLife committed
757
758


759
		//11. Restore parallel connectivity, global ids
Kirill Terekhov's avatar
Kirill Terekhov committed
760
		m->ResolveModification();
Kirill Terekhov's avatar
Kirill Terekhov committed
761
		//m->SynchronizeMarker(m->NewMarker(),CELL|FACE|EDGE|NODE,SYNC_BIT_OR);
762
        //ExchangeGhost(3,NODE); // Construct Ghost cells in 2 layers connected via nodes
763
764
		//12. Let the user update their data
		//todo: call back function for New() cells
Kirill Terekhov's avatar
Kirill Terekhov committed
765
		if( model ) model->Adaptation(*m);
766
		//13. Delete old elements of the mesh
Kirill Terekhov's avatar
Kirill Terekhov committed
767
		m->ApplyModification();
Kirill Terekhov's avatar
Kirill Terekhov committed
768
		
Kirill Terekhov's avatar
Kirill Terekhov committed
769
		//m->ExchangeGhost(1,NODE,m->NewMarker());
770
		//14. Done
Kirill Terekhov's avatar
Kirill Terekhov committed
771
        //cout << rank << ": Before end " << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
772
		m->EndModification();
773
		EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
774
		
775
776
777
		
		m->Save("after_refine"+std::to_string(fi)+".pvtk");
		std::cout << "Save after_refine"+std::to_string(fi)+".pvtk" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
778
		fi++;
Kirill Terekhov's avatar
Kirill Terekhov committed
779
		
780
		//ExchangeData(hanging_nodes,CELL | FACE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
781
        //m->ResolveSets();
782

Kirill Terekhov's avatar
Kirill Terekhov committed
783
784
785
786
787
788
        //m->BeginModification();
        //    m->ExchangeGhost(1,NODE,marker_new); // Construct Ghost cells in 2 layers connected via nodes
        //    m->ReleaseMarker(marker_new,CELL|FACE|EDGE|NODE);
		//m->ApplyModification();
    	//m->EndModification();
    	//PrintSet();
Kirill Terekhov's avatar
Kirill Terekhov committed
789
    	//m->ResolveSets();
Kirill Terekhov's avatar
Kirill Terekhov committed
790
		
Kirill Terekhov's avatar
Kirill Terekhov committed
791
		//restore face orientation
Kirill Terekhov's avatar
Kirill Terekhov committed
792
		//BUG: bad orientation not fixed automatically
793
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
794
795
796
797
798
799
800
801
		int nfixed = 0;
		for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) 
			if( !it->CheckNormalOrientation() )
			{
				it->FixNormalOrientation();
				nfixed++;
			}
			//std::cout << "Face " << it->LocalID() << " oriented incorrectly " << std::endl;
802
803
		if( nfixed ) REPORT_STR(rank << " fixed " << nfixed << " faces");
		EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
804
805
806
807
808
809
810
811
812
		/*
		for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
		{
			if( parent_set[*it] == InvalidHandle() )
				std::cout << m->GetProcessorRank() << " parent set not valid on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << std::endl;
			else if( GetHandleElementType(parent_set[*it]) != ESET )
				std::cout << m->GetProcessorRank() << " parent set is something else " << ElementTypeName(GetHandleElementType(parent_set[*it])) << ":" << GetHandleID(parent_set[*it]) << " on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << std::endl;
		}
		*/
Kirill Terekhov's avatar
Kirill Terekhov committed
813
814
		//std::cout << rank << " check parent_set at end" << std::endl;
		//CheckParentSet();
815
816

        //ExchangeData(hanging_nodes,CELL | FACE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
817
        //cout << rank << ": After end " << std::endl;
818
		//reorder element's data to free up space
819
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
820
		m->ReorderEmpty(CELL|FACE|EDGE|NODE);
821
		EXIT_BLOCK();
822
823
		//return number of refined cells
		call_counter++;
824
        EXIT_FUNC();
825
826
		return ret != 0;
	}
827
828
829
830
831
832
833
834
835
836
837
838

    void OperationMin(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
    {
        int value  =  (int)*((int*)data);


        if (value < element->Integer(tag))
        {
            element->Integer(tag) = value;
        }
        (void)size;
    }
839
	
840
841
    void AdaptiveMesh::SynchronizeIndicated(TagInteger& indicator)
    {
842
        if (m->GetProcessorsNumber() == 1) return;
843
        ENTER_FUNC();
844
        int rank = m->GetProcessorRank();
845
846

        // Check all sets. All elements in sets must be indicated. At first we check indicator in local processor, and second integrate data
847
848
        TagInteger tag_indicated = m->CreateTag("INDICATED",DATA_INTEGER,ESET,NONE,1);
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
849
        {
850
            ElementSet set = ElementSet(m,*it);
851
852
853
854
855
856
857
            set.Integer(tag_indicated) = 1;
            ElementSet::iterator p = set.Begin();
            while(p != set.End())
            {
                if (indicator[*p] == 0)
                {
                    tag_indicated[set] = 0;
858
                    //std::cout << rank << ": Set " << set.GetName() << " not all indicated" << endl;
859
860
861
862
863
864
865
                    break;
                }

                p++;
            }
        }

866
867
        m->ReduceData(tag_indicated,ESET,0,OperationMin);
        m->ExchangeData(tag_indicated,ESET,0);
868

869
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
870
871
872
873
874
875
876
877
878
            if (it->Integer(tag_indicated) == 0)
            {
                ElementSet::iterator p = it->Begin();
                while(p != it->End())
                {
                    p->Integer(indicator) = 0;
                    p++;
                }
            }
879
        EXIT_FUNC();
880
881
    }

882
883
	bool AdaptiveMesh::Coarse(TagInteger & indicator)
	{
884
885
		static int fi = 0;
        ENTER_FUNC();
Kirill Terekhov's avatar
Kirill Terekhov committed
886
        //return false;
887
888
889
890
		static int call_counter = 0;
		//return number of coarsened cells
		int ret = 0;
		//initialize tree structure
891
		ENTER_BLOCK();
892
		PrepareSet();
893
894
895
896
		EXIT_BLOCK();

		m->Save("before_coarse"+std::to_string(fi)+".pvtk");
		std::cout << "Save before_coarse"+std::to_string(fi)+".pvtk" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
897
		
898
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
899
900
901
902
903
904
905
906
907
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(parent_set,CELL,0);
		m->CheckCentroids(__FILE__,__LINE__);
		m->ResolveSets();
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(parent_set,CELL,0);
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
		m->CheckCentroids(__FILE__,__LINE__);
908
909
910
911
912
		CheckParentSet();
		EXIT_BLOCK();

		m->Save("before_coarse_parent"+std::to_string(fi)+".pvtk");
		std::cout << "Save before_coarse_parent"+std::to_string(fi)+".pvtk" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
913
		
914
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
915
		SynchronizeIndicated(indicator);
916
917
		EXIT_BLOCK();
		
918
919
		int schedule_counter = 1; //indicates order in which refinement will be scheduled
		int scheduled = 1, unscheduled = 0; //indicates that at least one element was scheduled on current sweep
920
921

		ENTER_BLOCK();
922
		//TagInteger coarsened = CreateTag("COARSENED",DATA_INTEGER,CELL,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
923
		TagInteger coarse_indicator = m->CreateTag("COARSE_INDICATOR",DATA_INTEGER,EDGE,NONE,1); //used to find on fine cells indicator on coarse cells
924
		//0. Extend indicator for sets, edges and faces
Kirill Terekhov's avatar
Kirill Terekhov committed
925
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
926
927
928
929
930
931
932
933
		while(scheduled || unscheduled)
		{
			// rules
			// a) If there is adjacent finer edge that is not marked for coarsening
			// then this cell should not be coarsened
			// b) If there is adjacent coarser cell, then this cell should be coarsened
			// first
			//0.Communicate indicator - it may be not synced
Kirill Terekhov's avatar
Kirill Terekhov committed
934
			m->ExchangeData(indicator,CELL,0);
935
936
937
938
939
940
941
			//1. Mark each adjacent face/edge for coarsement schedule
			// problem: should mark so that if every adjacent cell is coarsened
			// then adjacent face/edge are also coarsened
			for(ElementType etype = EDGE; etype <= FACE; etype = NextElementType(etype))
			{
				//for(Storage::integer it = 0; it < LastLocalID(etype); ++it) if( isValidElement(etype,it) )
				//	indicator[ElementByLocalID(etype,it)] = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
942
				for(Storage::integer it = 0; it < m->LastLocalID(etype); ++it) if( m->isValidElement(etype,it) )
943
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
944
					Element e = m->ElementByLocalID(etype,it);
945
946
947
948
					ElementArray<Cell> adj = e.getCells();
					indicator[e] = INT_MAX;
					for(ElementArray<Element>::size_type kt = 0; kt < adj.size(); ++kt)
						if( level[e] == level[adj[kt]]) indicator[e] = std::min(indicator[e],indicator[adj[kt]]);
949
					//if (indicator[e] == INT_MAX) cout << rank << ": " << ElementTypeName(e.GetElementType()) << e.GlobalID() << endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
950
					//assert(indicator[e] != INT_MAX);
951
952
953
				}
			}
			//2.Communicate indicator on faces and edges
Kirill Terekhov's avatar
Kirill Terekhov committed
954
			m->ReduceData(indicator,FACE|EDGE,0,ReduceMin);
955
			m->ExchangeData(indicator,FACE|EDGE,0);
956
957
958
			//3.If there is adjacent finer edge that are not marked for coarsening
			// then this cell should not be coarsened
			unscheduled = scheduled = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
959
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
960
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
961
				Cell c = m->CellByLocalID(it);
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
				if( indicator[c] )
				{
					ElementArray<Edge> edges = c.getEdges();
					for(ElementArray<Edge>::size_type kt = 0; kt < edges.size(); ++kt)
					{
						if( level[edges[kt]] > level[c] && indicator[edges[kt]] == 0 )
						{
							indicator[c] = 0;
							unscheduled++;
						}
					}
				}
			}
			//4. Propogate coarsement info over set tree to detect valid coarsenings.
			// go down over sets, if set does not have children and all of the cells
			// of the set are marked for coarsening, then mark the set for coarsement
			// otherwise unmark.
			// Unmark all cells that are not to be coarsened
Kirill Terekhov's avatar
Kirill Terekhov committed
980
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
981
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
982
				Cell c = m->CellByLocalID(it);
983
                //if (!isValidElement(c.GetHandle())) continue;
984
985
				if( indicator[c] )
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
986
					ElementSet parent(m,parent_set[c]);
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
					//intermediate cell may not be coarsened
					//root set may not have coarsening cells
					if( parent.HaveChild() || !parent.HaveParent() )
					{
						indicator[c] = 0;
						unscheduled++;
					}
					else
					{
						Storage::integer schedule_first = 0;
						bool check = true;
						//check that all elements of the set are to be coarsened
						for(ElementSet::iterator it = parent.Begin(); it != parent.End(); ++it)
						{
							check &= (indicator[it->self()] != 0);
							schedule_first = std::max(schedule_first,indicator[it->self()]);
						}
						if(!check)
						{
							indicator[c] = 0;
							unscheduled++;
						}
						else if( indicator[c] != schedule_first )
						{
							indicator[c] = schedule_first;
							unscheduled++;
						}
					}
				}
			}
			//5.If there is an adjacent coarser element to be refined, then
			//   this one should be scheduled to be refined first
			//a) clean up coarse indicator tag
Kirill Terekhov's avatar
Kirill Terekhov committed
1020
1021
			for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
				coarse_indicator[m->EdgeByLocalID(it)] = 0;
1022
			//b) each cell mark it's finer edges with cell's schedule
Kirill Terekhov's avatar
Kirill Terekhov committed
1023
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
1024
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1025
				Cell c = m->CellByLocalID(it);
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
				if( indicator[c] )
				{
					ElementArray<Element> adj = c.getAdjElements(EDGE);
					for(ElementArray<Element>::size_type kt = 0; kt < adj.size(); ++kt)
					{
						if( level[adj[kt]] > level[c] ) //only finer edges
							coarse_indicator[adj[kt]] = std::max(coarse_indicator[adj[kt]],indicator[c]);
					}
				}
			}
			//c) data reduction to get maximum over mesh partition
Kirill Terekhov's avatar
Kirill Terekhov committed
1037
			m->ReduceData(coarse_indicator,EDGE,0,ReduceMax);
1038
			m->ExchangeData(coarse_indicator,EDGE,0);
1039
			//d) look from cells if any edge is coarsened earlier
Kirill Terekhov's avatar
Kirill Terekhov committed
1040
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
1041
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1042
				Cell c = m->CellByLocalID(it);
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
				if( indicator[c] )
				{
					ElementArray<Element> adj = c.getAdjElements(EDGE);
					for(ElementArray<Element>::size_type kt = 0; kt < adj.size(); ++kt)
					{
						if( level[c] == level[adj[kt]] && //do not look from coarse cell onto finer edge
						    indicator[c] <= coarse_indicator[adj[kt]])
						{
							indicator[c] = coarse_indicator[adj[kt]]+1;
							scheduled++;
						}
					}
				}
			}
1057
			m->ExchangeData(indicator,CELL|FACE|EDGE,0);
1058
			//5.Go back to 1 until no new elements scheduled
Kirill Terekhov's avatar
Kirill Terekhov committed
1059
1060
			scheduled = m->Integrate(scheduled);
			unscheduled = m->Integrate(unscheduled);
1061
1062
1063
			if( scheduled ) schedule_counter++;
		}
		//cleanup
Kirill Terekhov's avatar
Kirill Terekhov committed
1064
		coarse_indicator = m->DeleteTag(coarse_indicator);
1065
		EXIT_BLOCK();
1066
		//Make schedule which elements should be refined earlier.
1067
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
1068
		m->BeginModification();
1069
1070
1071
1072
1073
1074
		while(schedule_counter)
		{
			//unite cells
			//should find and set hanging nodes on faces
			//find single node at the center, all other nodes,
			//adjacent over edge are hanging nodes
Kirill Terekhov's avatar
Kirill Terekhov committed
1075
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
1076
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1077
				Cell c = m->CellByLocalID(it);
1078
1079
1080
				if( !c.Hidden() && indicator[c] == schedule_counter )
				{
					//this set contains all the cells to be united
Kirill Terekhov's avatar
Kirill Terekhov committed
1081
1082
					ElementSet parent(m,parent_set[c]);
					ElementArray<Cell> unite_cells(m,parent.Size());
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
					//unmark indicator to prevent coarsement with next element
					Storage::integer kt = 0;
					for(ElementSet::iterator jt = parent.Begin(); jt != parent.End(); ++jt)
					{
						unite_cells[kt++] = jt->getAsCell();
						indicator[jt->self()] = 0; //prevent algorithm from visiting again
					}
					//find a node common to all the cells
					ElementArray<Node> center_node = unite_cells[0].getNodes();
					for(kt = 1; kt < unite_cells.size(); ++kt)
						center_node.Intersect(unite_cells[kt].getNodes());
					//only one should be found
					assert(center_node.size() == 1);
					ElementArray<Node> hanging =