amesh.cpp 57 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


#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();}
11
#define ENTER_BLOCK() { double btime = Timer(); m->WriteTab(m->out_time) << "<FUNCTION name=\"" << __FUNCTION__ << ":" << NameSlash(__FILE__) << ":" << __LINE__ << "\" id=\"func" << m->GetFuncID()++ << "\">" << std::endl; m->Enter();
12
13
14
15
16
17
18
#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) {}
19
20
#define ENTER_BLOCK()
#define EXIT_BLOCK()
21
22
23
24
25
#define ENTER_FUNC() {}
#define EXIT_FUNC() {}
#define EXIT_FUNC_DIE()  {}
#endif

26
27
28
29
30
31
32
33
static std::string NameSlash(std::string input)
{
	for(unsigned l = input.size(); l > 0; --l)
		if( input[l-1] == '/' || input[l-1] == '\\' )
			return std::string(input.c_str() + l);
	return input;
}

34
35

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

Kirill Terekhov's avatar
Kirill Terekhov committed
38
39
//from inmost
//std::string base64_encode(unsigned char const* buf, unsigned int bufLen);
40

Kirill Terekhov's avatar
Kirill Terekhov committed
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
77
78
79
80
81
82
83
84
85
86
	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;
	}
87

88
89
90
/// todo:
/// 1. coarsment
/// 2. strategy for faces/edges with faults
Kirill Terekhov's avatar
Kirill Terekhov committed
91
/// 3. geom model supportbn
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
/// 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
119

120
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
121
    void AdaptiveMesh::PrintSetLocal(std::string offset, ElementSet it, std::stringstream& ss)
SilverLife's avatar
SilverLife committed
122
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
123
124
        std::stringstream ss1;
        ss1 << offset << rank << ": Set : " << std::setw(5) << it.GetName() << " ";
SilverLife's avatar
SilverLife committed
125
126
        for (int i = ss1.str().length(); i < 23; i++) ss1 << " ";
        ss << ss1.str();
Kirill Terekhov's avatar
Kirill Terekhov committed
127
        ss << std::setw(6);
SilverLife's avatar
SilverLife committed
128
129
130
131
132
        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";

133
        ss << " tag_owner (" << it.IntegerDF(m->OwnerTag()) << ")";
SilverLife's avatar
SilverLife committed
134
135
136

        //ss << "   level (" << level[it.self()] << ")  ";
        ss << " tag_processors (";
Kirill Terekhov's avatar
Kirill Terekhov committed
137
        std::stringstream ss2;
138
        Storage::integer_array arr = it.IntegerArrayDV(m->ProcessorsTag());
SilverLife's avatar
SilverLife committed
139
140
        for (int i = 0; i < arr.size(); i++)
            ss2 << arr[i] << " ";
Kirill Terekhov's avatar
Kirill Terekhov committed
141
        ss << std::setw(5) << ss2.str() <<")";
SilverLife's avatar
SilverLife committed
142
143
144
145
146
147
148
149
        

        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
150
            std::string type = "unknw";
SilverLife's avatar
SilverLife committed
151
152
153
154
            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
155
            ss << type << "-" << std::setw(2) << p->GlobalID() << " ";
SilverLife's avatar
SilverLife committed
156
157
            p++;
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
158
        ss << std::endl;
SilverLife's avatar
SilverLife committed
159
160
161
162
163
164

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

167
	/*
SilverLife's avatar
SilverLife committed
168
169
    void AdaptiveMesh::PrintSet()
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
170
        std::stringstream ss;
171
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
SilverLife's avatar
SilverLife committed
172
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
173
            //if (it->HaveParent()) continue;
174
            PrintSetLocal("",ElementSet(m,*it),ss);
SilverLife's avatar
SilverLife committed
175
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
176
        std::cout << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
177
    }
178
	*/
SilverLife's avatar
SilverLife committed
179

180
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
181
    void PrintRefs(std::ostream& os, Storage::reference_array refs)
182
183
184
    {
        for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
185
            std::string type = "unknw";
186
187
188
189
190
191
192
            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
193
            os << "(" << xyz[0] << "," << xyz[1] << "," << xyz[2] <<")" <<  std::endl;
194
195
        }
    }
196
	*/
197

198
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
199
    void AdaptiveMesh::PrintMesh(std::ostream& os, int cell, int face, int edge, int node)
SilverLife's avatar
SilverLife committed
200
201
    {
        if (cell + face + edge + node == 0) return;
Kirill Terekhov's avatar
Kirill Terekhov committed
202
203
        std::stringstream ss;
        ss << "================= " << rank << " =====================" << std::endl;
SilverLife's avatar
SilverLife committed
204
205
        if (cell)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
206
            ss << "Cells: " << m->NumberOfCells() <<  std::endl;
207
            for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) 
SilverLife's avatar
SilverLife committed
208
            {
209
                ss << rank << ": " << it->GlobalID() << " - " << it->LocalID() << " - ";            
SilverLife's avatar
SilverLife committed
210
211
212
                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
213
                ss << std::endl;
SilverLife's avatar
SilverLife committed
214
215
216
217
218
            }
        }

        if (face)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
219
            ss << "Faces: " << m->NumberOfFaces() <<  std::endl;
220
            for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) 
SilverLife's avatar
SilverLife committed
221
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
222
223
                ss << rank << ": " << std::setw(2) << it->LocalID() << " " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
224
225
226
227
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";

228
229
230

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

233
                ss << "  " << m->GetMarker(*it,m->NewMarker());
234
235
236
237

                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
238
                    ss << std::setw(2) << node->GlobalID() << " ";
239
                ss << ")";
Kirill Terekhov's avatar
Kirill Terekhov committed
240
                ss << std::endl;
SilverLife's avatar
SilverLife committed
241
242
243
244
245
            }
        }

        if (edge)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
246
            ss << "Edges: " << m->NumberOfEdges() <<  std::endl;
247
            for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it) 
SilverLife's avatar
SilverLife committed
248
249
250
251
252
            {
                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
253
                ss << std::endl;
SilverLife's avatar
SilverLife committed
254
255
256
257
258
            }
        }

        if (node)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
259
            ss << "Nodes:" << std::endl;
260
            for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) 
SilverLife's avatar
SilverLife committed
261
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
262
263
                ss << rank << ": " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
264
265
266
267
268
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";

                {
269
                    ss << "  " << m->GetMarker(*it,m->NewMarker());
270
271

                    ss << "(" << 
Kirill Terekhov's avatar
Kirill Terekhov committed
272
273
274
                        std::setw(3) << it->RealArray(m->CoordsTag())[0] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[1] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[2] << ")";
275

SilverLife's avatar
SilverLife committed
276
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
277
                ss << std::endl;
SilverLife's avatar
SilverLife committed
278
279
280
            }
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
281
282
        ss << "=========================================" << std::endl;
        os << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
283
    }
284
	*/
SilverLife's avatar
SilverLife committed
285

286
	/*
SilverLife's avatar
SilverLife committed
287
288
289
290
291
    void AdaptiveMesh::UpdateStatus()
    {

        for(ElementType mask = CELL; mask >= NODE; mask = PrevElementType(mask))
        {
292
            for(Mesh::iteratorElement it = m->BeginElement(mask); it != m->EndElement(); it++)
SilverLife's avatar
SilverLife committed
293
294
295
296
297
298
299
300
301
            {
                int stat = 0;
                if (it->GetStatus() == Element::Shared) stat = 1;
                else if (it->GetStatus() == Element::Ghost)  stat = 2;

                tag_status[it->self()] = stat;
            }
        }
    }
302
	*/
SilverLife's avatar
SilverLife committed
303
    
304
305
	
	void AdaptiveMesh::PrintSet(std::ostream & fout, ElementSet set)
SilverLife's avatar
SilverLife committed
306
    {
307
308
309
		fout << "set " << set.GetName();
		fout << " size " << set.Size();
		if(set.HaveParent()) fout << " parent " << set.GetParent().GetName();
310
311
312
		fout << " elements ";
		for(ElementSet::iterator it = set.Begin(); it != set.End(); ++it)
			fout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << ":" << it->GlobalID() << " ";
313
		fout << " children ";
SilverLife's avatar
SilverLife committed
314
		for(ElementSet child = set.GetChild(); child.isValid(); child = child.GetSibling())
315
316
            fout << child.GetName() << " ";
		fout << std::endl;
SilverLife's avatar
SilverLife committed
317
    }
318
	
SilverLife's avatar
SilverLife committed
319

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

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

438
439
		//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
440
		
441
442
443
444
445
		//ENTER_BLOCK();
		//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
		//	if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
		//EXIT_BLOCK();
		
446
		ENTER_BLOCK();
447
		m->ResolveSets();
448
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
449
		m->ExchangeData(parent_set,CELL,0);
450
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
451
		m->ResolveSets();
452
		//m->CheckCentroids(__FILE__,__LINE__);
453
		//m->ExchangeData(parent_set,CELL,0);
454
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
455
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
456
457
		//m->CheckCentroids(__FILE__,__LINE__);
		//CheckParentSet();
458
		EXIT_BLOCK();
459
460
461
462
463
		
		//ENTER_BLOCK();
		//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
		//	if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
		//EXIT_BLOCK();
464

465
466
		//m->Save("before_refine_parent"+std::to_string(fi)+".pvtk");
		//std::cout << "Save before_refine_parent"+std::to_string(fi)+".pvtk" << std::endl;
467
468
469
		

		
470
471
		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
472
		ENTER_BLOCK();
473
		//0. Extend indicator for edges and faces
Kirill Terekhov's avatar
Kirill Terekhov committed
474
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
475
476
		while(scheduled)
		{
477
			REPORT_VAL("scheduled",scheduled);
478
			//1.Communicate indicator - it may be not synced
Kirill Terekhov's avatar
Kirill Terekhov committed
479
			m->ExchangeData(indicator,CELL,0);
480
481
			//2.Propogate indicator down to the faces,edges
			//  select schedule for them
482
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
483
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
484
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
485
				Cell c = m->CellByLocalID(it);
486
487
488
489
490
491
492
493
494
495
				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
					}
				}
			}
496
			EXIT_BLOCK();
497
			//3.Communicate indicator on faces and edges
Kirill Terekhov's avatar
Kirill Terekhov committed
498
			m->ExchangeData(indicator,FACE|EDGE,0);
499
500
501
502
			//4.Check for each cell if there is
			//  any hanging node with adjacent in a need to refine,
			//  schedule for refinement earlier.
			scheduled = 0;
503
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
504
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
505
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
506
				Cell c = m->CellByLocalID(it);
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
				//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;
						}
					}
				}
			}
528
			EXIT_BLOCK();
529
			//5.Go back to 1 until no new elements scheduled
Kirill Terekhov's avatar
Kirill Terekhov committed
530
			scheduled = m->Integrate(scheduled);
531
532
			if( scheduled ) schedule_counter++;
		}
533
		m->ExchangeData(indicator,CELL | FACE | EDGE,0);
534
		EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
535
		//m->Save("indicator.pmf");
536
537
538
539
540
541
542
543
		
		//ENTER_BLOCK();
		//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
		//	if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
		//EXIT_BLOCK();
		
		//if( !Element::CheckConnectivity(m) ) std::cout << __FILE__ << ":" << __LINE__ << " broken connectivity" << std::endl;
		//m->CheckCentroids(__FILE__,__LINE__);
544
		//6.Refine
545
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
546
		m->BeginModification();
547
548
549
550
		while(schedule_counter)
		{
			Storage::real xyz[3] = {0,0,0};
			//7.split all edges of the current schedule
551
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
552
			for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
553
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
554
				Edge e = m->EdgeByLocalID(it);
555
556
				if( !e.Hidden() && indicator[e] == schedule_counter )
				{
557
558
559
560
561
562
563
					//ENTER_BLOCK();
					//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
					//	if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
					//EXIT_BLOCK();
					//REPORT_STR("split edge " << e.LocalID() << " nodes " << e.getBeg().LocalID() << "," << e.getEnd().LocalID() << " level " << level[e] << " lc size " << m->LowConn(e.GetHandle()).size() );
					//ElementArray<Node> nodes = e.getNodes();
					//for(int q = 0; q < nodes.size(); ++q) REPORT_STR("node " << nodes[q].GetHandle() << " " << nodes[q].LocalID() << (nodes[q].Hidden()?" hidden " : " good ") );
564
565
566
					//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
567
					for(Storage::integer d = 0; d < m->GetDimensions(); ++d)
568
569
570
						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
571
					Node n = m->CreateNode(xyz);
572
					//set increased level for new node
573
					level[n] = level[e.getBeg()] = level[e.getEnd()] = level[e]+1;
574
575
576
577
					//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
578
					ElementArray<Edge> new_edges = Edge::SplitEdge(e,ElementArray<Node>(m,1,n.GetHandle()),0);
579
580
					//set increased level for new edges
					level[new_edges[0]] = level[new_edges[1]] = level[e]+1;
581
582
583
584
585
586
587
					
					//for(int q = 0; q < 2; ++q)
					//{
					//	REPORT_STR("new edges["<<q<<"]" << new_edges[q].LocalID() << " nodes " << new_edges[q].getBeg().LocalID() << "," << new_edges[q].getEnd().LocalID() << " level " << level[new_edges[q]]);
					//}
					
					//if( !Element::CheckConnectivity(m) ) std::cout << __FILE__ << ":" << __LINE__ << " broken connectivity" << std::endl;
588
589
				}
			}
590
			EXIT_BLOCK();
591
			//8.split all faces of the current schedule, using hanging nodes on edges
592
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
593
			for(Storage::integer it = 0; it < m->FaceLastLocalID(); ++it) if( m->isValidFace(it) )
594
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
595
				Face f = m->FaceByLocalID(it);
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
				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
611
					Node n = m->CreateNode(xyz);
612
					//set increased level for the new node
613
					level[n] = level[f]+1;
614
615
616
					//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
617
618
					ElementArray<Node> edge_nodes(m,2); //to create new edges
					ElementArray<Edge> hanging_edges(m,face_hanging_nodes.size());
619
620
621
622
					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
623
						hanging_edges[kt] = m->CreateEdge(edge_nodes).first;
624
625
626
627
628
629
630
631
632
633
						//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;
				}
			}
634
			EXIT_BLOCK();
635
			//this tag helps recreate internal face
Kirill Terekhov's avatar
Kirill Terekhov committed
636
			TagReferenceArray internal_face_edges = m->CreateTag("INTERNAL_FACE_EDGES",DATA_REFERENCE,NODE,NODE,4);
637
			//this marker helps detect edges of current cell only
Kirill Terekhov's avatar
Kirill Terekhov committed
638
			MarkerType mark_cell_edges = m->CreateMarker();
639
			//this marker helps detect nodes hanging on edges of unrefined cell
Kirill Terekhov's avatar
Kirill Terekhov committed
640
			MarkerType mark_hanging_nodes = m->CreateMarker();
641
			//9.split all cells of the current schedule
642
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
643
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
644
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
645
				Cell c = m->CellByLocalID(it);
646
647
648
649
650
651
652
653
654
655
656
				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
657
					Node n = m->CreateNode(xyz);
658
					//set increased level for the new node
659
					level[n] = level[c]+1;
660
661
662
663
664
					//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
665
666
					ElementArray<Node> edge_nodes(m,2);
					ElementArray<Edge> edges_to_faces(m,cell_hanging_nodes.size());
667
668
669
670
671
672
					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
673
						edges_to_faces[kt] = m->CreateEdge(edge_nodes).first;
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
						//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
709
					ElementArray<Face> internal_faces(m,edge_hanging_nodes.size());
710
711
712
713
714
715
716
717
718
719
					//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
720
						internal_faces[kt] = m->CreateFace(ElementArray<Edge>(m,face_edges.begin(),face_edges.end())).first;
721
722
723
724
725
						//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
726
727
728
729
730
731
732
					//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;
					//	
					//}
733
734
					//split the cell
					ElementArray<Cell> new_cells = Cell::SplitCell(c,internal_faces,0);
735
					std::sort(new_cells.begin(),new_cells.end(),Mesh::CentroidComparator(m));
736
					//retrive parent set
Kirill Terekhov's avatar
Kirill Terekhov committed
737
					ElementSet parent(m,parent_set[c]);
738
					//create set corresponding to old coarse cell
Kirill Terekhov's avatar
Kirill Terekhov committed
739
740
					Storage::real cnt[3];
					c.Centroid(cnt);
741
					std::stringstream set_name;
Kirill Terekhov's avatar
Kirill Terekhov committed
742
					//set_name << parent.GetName() << "_C" << c.GlobalID(); //rand may be unsafe
Kirill Terekhov's avatar
Kirill Terekhov committed
743
					if( parent == root )
744
						set_name << "AM_R" << set_id[c];
Kirill Terekhov's avatar
Kirill Terekhov committed
745
					else
746
						set_name << parent.GetName() << "C" << set_id[c];
747
					//set_name << base64_encode_((unsigned char *)cnt,3*sizeof(double)/sizeof(unsigned char));
Kirill Terekhov's avatar
Kirill Terekhov committed
748
					
Kirill Terekhov's avatar
Kirill Terekhov committed
749
750
751
752
753
754
755
756
757
758
					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
759
						exit(-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
760
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
761
					
Kirill Terekhov's avatar
Kirill Terekhov committed
762
					ElementSet cell_set = m->CreateSetUnique(set_name.str()).first;
763
					//cell_set->SetExchange(ElementSet::SYNC_ELEMENTS_ALL);
764
765
766
767
					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
768
						set_id[new_cells[kt]] = kt;
769
770
771
772
						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
773
774
775
776
777
778
779
780
781
782
					/*
					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() )
783
					parent.AddChild(cell_set);
Kirill Terekhov's avatar
Kirill Terekhov committed
784
					//else assert(cell_set->GetParent() == parent);
785
786
787
788
					//increment number of refined cells
					ret++;
				}
			}
789
			EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
790
791
792
			m->ReleaseMarker(mark_hanging_nodes);
			m->ReleaseMarker(mark_cell_edges);
			m->DeleteTag(internal_face_edges);
793
794
795
			//10.jump to later schedule, and go to 7.
			schedule_counter--;
		}
SilverLife's avatar
SilverLife committed
796

797
		//free created tag
Kirill Terekhov's avatar
Kirill Terekhov committed
798
		m->DeleteTag(indicator,FACE|EDGE);
SilverLife's avatar
SilverLife committed
799
800


801
		//11. Restore parallel connectivity, global ids
Kirill Terekhov's avatar
Kirill Terekhov committed
802
		m->ResolveModification();
Kirill Terekhov's avatar
Kirill Terekhov committed
803
		//m->SynchronizeMarker(m->NewMarker(),CELL|FACE|EDGE|NODE,SYNC_BIT_OR);
804
        //ExchangeGhost(3,NODE); // Construct Ghost cells in 2 layers connected via nodes
805
806
		//12. Let the user update their data
		//todo: call back function for New() cells
Kirill Terekhov's avatar
Kirill Terekhov committed
807
		if( model ) model->Adaptation(*m);
808
		//13. Delete old elements of the mesh
Kirill Terekhov's avatar
Kirill Terekhov committed
809
		m->ApplyModification();
Kirill Terekhov's avatar
Kirill Terekhov committed
810
		
Kirill Terekhov's avatar
Kirill Terekhov committed
811
		//m->ExchangeGhost(1,NODE,m->NewMarker());
812
		//14. Done
Kirill Terekhov's avatar
Kirill Terekhov committed
813
        //cout << rank << ": Before end " << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
814
		m->EndModification();
815
		EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
816
		
817
		
818
819
		//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
820
		fi++;
Kirill Terekhov's avatar
Kirill Terekhov committed
821
		
822
		//ExchangeData(hanging_nodes,CELL | FACE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
823
        //m->ResolveSets();
824

Kirill Terekhov's avatar
Kirill Terekhov committed
825
826
827
828
829
830
        //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
831
    	//m->ResolveSets();
Kirill Terekhov's avatar
Kirill Terekhov committed
832
		
Kirill Terekhov's avatar
Kirill Terekhov committed
833
		//restore face orientation
Kirill Terekhov's avatar
Kirill Terekhov committed
834
		//BUG: bad orientation not fixed automatically
835
836
		
		/*
837
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
838
839
840
841
842
843
844
845
		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;
846
847
		if( nfixed ) REPORT_STR(rank << " fixed " << nfixed << " faces");
		EXIT_BLOCK();
848
849
		 */
		
Kirill Terekhov's avatar
Kirill Terekhov committed
850
851
852
853
854
855
856
857
858
		/*
		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
859
860
		//std::cout << rank << " check parent_set at end" << std::endl;
		//CheckParentSet();
861
862

        //ExchangeData(hanging_nodes,CELL | FACE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
863
        //cout << rank << ": After end " << std::endl;
864
		//reorder element's data to free up space
865
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
866
		m->ReorderEmpty(CELL|FACE|EDGE|NODE);
867
		EXIT_BLOCK();
868
869
		//return number of refined cells
		call_counter++;
870
        EXIT_FUNC();
871
872
		return ret != 0;
	}
873
874
875
876
877
878
879
880
881
882
883
884

    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;
    }
885
	
886
887
    void AdaptiveMesh::SynchronizeIndicated(TagInteger& indicator)
    {
888
        if (m->GetProcessorsNumber() == 1) return;
889
        ENTER_FUNC();
890
        int rank = m->GetProcessorRank();
891
892

        // Check all sets. All elements in sets must be indicated. At first we check indicator in local processor, and second integrate data
893
        TagInteger tag_indicated = m->CreateTag("INDICATED",DATA_INTEGER,ESET,NONE,1);
894
		ENTER_BLOCK();
895
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
896
        {
897
            ElementSet set = ElementSet(m,*it);
898
899
900
901
902
903
904
            set.Integer(tag_indicated) = 1;
            ElementSet::iterator p = set.Begin();
            while(p != set.End())
            {
                if (indicator[*p] == 0)
                {
                    tag_indicated[set] = 0;
905
                    //std::cout << rank << ": Set " << set.GetName() << " not all indicated" << endl;
906
907
908
909
910
911
                    break;
                }

                p++;
            }
        }
912
		EXIT_BLOCK();
913
914
        m->ReduceData(tag_indicated,ESET,0,OperationMin);
        m->ExchangeData(tag_indicated,ESET,0);
915
		ENTER_BLOCK();
916
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
917
918
919
920
921
922
923
924
925
            if (it->Integer(tag_indicated) == 0)
            {
                ElementSet::iterator p = it->Begin();
                while(p != it->End())
                {
                    p->Integer(indicator) = 0;
                    p++;
                }
            }
926
		EXIT_BLOCK();
927
        EXIT_FUNC();
928
929
    }

930
931
	bool AdaptiveMesh::Coarse(TagInteger & indicator)
	{
932
		std::string file;
933
		ENTER_FUNC();
Kirill Terekhov's avatar
Kirill Terekhov committed
934
        //return false;
935
936
937
938
		static int call_counter = 0;
		//return number of coarsened cells
		int ret = 0;
		//initialize tree structure
939
		ENTER_BLOCK();
940
		PrepareSet();
941
942
		EXIT_BLOCK();

943
944
		//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
945
		
946
		ENTER_BLOCK();
947
		//m->CheckCentroids(__FILE__,__LINE__);
948
		m->ResolveSets();
Kirill Terekhov's avatar
Kirill Terekhov committed
949
		m->ExchangeData(parent_set,CELL,0);
950
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
951
		m->ResolveSets();
952
		//m->CheckCentroids(__FILE__,__LINE__);
953
		//m->ExchangeData(parent_set,CELL,0);
954
		//m->CheckCentroids(__FILE__,__LINE__);
955
		//m->ExchangeData(hanging_nodes,CELL | FACE,0);
956
957
		//m->CheckCentroids(__FILE__,__LINE__);
		//CheckParentSet();
958
959
		EXIT_BLOCK();

960
961
962
963
964
965
966
		//m->Save("before_coarse_parent"+std::to_string(fi)+".pvtk");
		//std::cout << "Save before_coarse_parent"+std::to_string(fi)+".pvtk" << std::endl;
		
		//ENTER_BLOCK();
		//SynchronizeIndicated(indicator);
		//EXIT_BLOCK();
		
Kirill Terekhov's avatar
Kirill Terekhov committed
967
		
968
		
969
970
		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
971
972

		ENTER_BLOCK();
973
		//TagInteger coarsened = CreateTag("COARSENED",DATA_INTEGER,CELL,NONE,1);
974
		TagInteger coarse_indicator = m->CreateTag("COARSE_INDICATOR",DATA_INTEGER,ESET|EDGE,ESET,1); //used to find on fine cells indicator on coarse cells
975
		//0. Extend indicator for sets, edges and faces
976
977
978
979
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,ESET|FACE|EDGE,ESET,1);
		std::vector<Tag> indicators;
		indicators.push_back(coarse_indicator);
		indicators.push_back(indicator);
980
981
982
983
		//int mi = 0;
		//file = "w"+std::to_string(mi++)+".pvtk";
		//m->Save(file);
		//std::cout << "Save " << file << " " << __FILE__ << ":" << __LINE__ << std::endl;
984
985
		while(scheduled || unscheduled)
		{
986
987
			REPORT_VAL("scheduled",scheduled);
			REPORT_VAL("unscheduled",unscheduled);
988
989
990
991
992
993
			// 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
994
			m->ExchangeData(indicator,CELL,0);
995
996
997
			//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
998
			ENTER_BLOCK();
999
1000
1001
1002
			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
1003
				for(Storage::integer it = 0; it < m->LastLocalID(etype); ++it) if( m->isValidElement(etype,it) )
1004
				{
Kirill Terekhov's avatar
Kirill Terekhov committed