amesh.cpp 48.8 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
19
20
21
22
23
#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

24
25
26
27
28
29
30
31
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;
}

32
33

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

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

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

86
87
88
/// todo:
/// 1. coarsment
/// 2. strategy for faces/edges with faults
Kirill Terekhov's avatar
Kirill Terekhov committed
89
/// 3. geom model supportbn
90
91
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
/// 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
117

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

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

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

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

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

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

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

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

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

226
227
228

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

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

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

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

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

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

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

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

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

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

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

                tag_status[it->self()] = stat;
            }
        }
    }
300
	*/
SilverLife's avatar
SilverLife committed
301
    
302
303
	
	void AdaptiveMesh::PrintSet(std::ostream & fout, ElementSet set)
SilverLife's avatar
SilverLife committed
304
    {
305
306
307
308
		fout << "set " << set.GetName();
		fout << " size " << set.Size();
		if(set.HaveParent()) fout << " parent " << set.GetParent().GetName();
		fout << " children ";
SilverLife's avatar
SilverLife committed
309
		for(ElementSet child = set.GetChild(); child.isValid(); child = child.GetSibling())
310
311
            fout << child.GetName() << " ";
		fout << std::endl;
SilverLife's avatar
SilverLife committed
312
    }
313
	
SilverLife's avatar
SilverLife committed
314

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

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

433
434
		//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
435
		
436
		ENTER_BLOCK();
437
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
438
		m->ExchangeData(parent_set,CELL,0);
439
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
440
		m->ResolveSets();
441
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
442
		m->ExchangeData(parent_set,CELL,0);
443
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
444
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
445
446
		//m->CheckCentroids(__FILE__,__LINE__);
		//CheckParentSet();
447
448
		EXIT_BLOCK();

449
450
		//m->Save("before_refine_parent"+std::to_string(fi)+".pvtk");
		//std::cout << "Save before_refine_parent"+std::to_string(fi)+".pvtk" << std::endl;
451
452
453
		

		
454
455
		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
456
		ENTER_BLOCK();
457
		//0. Extend indicator for edges and faces
Kirill Terekhov's avatar
Kirill Terekhov committed
458
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
459
460
		while(scheduled)
		{
461
			REPORT_VAL("scheduled",scheduled);
462
			//1.Communicate indicator - it may be not synced
Kirill Terekhov's avatar
Kirill Terekhov committed
463
			m->ExchangeData(indicator,CELL,0);
464
465
			//2.Propogate indicator down to the faces,edges
			//  select schedule for them
466
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
467
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
468
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
469
				Cell c = m->CellByLocalID(it);
470
471
472
473
474
475
476
477
478
479
				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
					}
				}
			}
480
			EXIT_BLOCK();
481
			//3.Communicate indicator on faces and edges
Kirill Terekhov's avatar
Kirill Terekhov committed
482
			m->ExchangeData(indicator,FACE|EDGE,0);
483
484
485
486
			//4.Check for each cell if there is
			//  any hanging node with adjacent in a need to refine,
			//  schedule for refinement earlier.
			scheduled = 0;
487
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
488
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
489
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
490
				Cell c = m->CellByLocalID(it);
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
				//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;
						}
					}
				}
			}
512
			EXIT_BLOCK();
513
			//5.Go back to 1 until no new elements scheduled
Kirill Terekhov's avatar
Kirill Terekhov committed
514
			scheduled = m->Integrate(scheduled);
515
516
			if( scheduled ) schedule_counter++;
		}
517
		m->ExchangeData(indicator,CELL | FACE | EDGE,0);
518
		EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
519
		//m->Save("indicator.pmf");
520
		//6.Refine
521
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
522
		m->BeginModification();
523
524
525
526
		while(schedule_counter)
		{
			Storage::real xyz[3] = {0,0,0};
			//7.split all edges of the current schedule
527
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
528
			for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
529
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
530
				Edge e = m->EdgeByLocalID(it);
531
532
533
534
535
				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
536
					for(Storage::integer d = 0; d < m->GetDimensions(); ++d)
537
538
539
						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
540
					Node n = m->CreateNode(xyz);
541
					//set increased level for new node
542
					level[n] = level[e.getBeg()] = level[e.getEnd()] = level[e]+1;
543
544
545
546
					//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
547
					ElementArray<Edge> new_edges = Edge::SplitEdge(e,ElementArray<Node>(m,1,n.GetHandle()),0);
548
549
550
551
					//set increased level for new edges
					level[new_edges[0]] = level[new_edges[1]] = level[e]+1;
				}
			}
552
			EXIT_BLOCK();
553
			//8.split all faces of the current schedule, using hanging nodes on edges
554
			ENTER_BLOCK();
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
						//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;
				}
			}
596
			EXIT_BLOCK();
597
			//this tag helps recreate internal face
Kirill Terekhov's avatar
Kirill Terekhov committed
598
			TagReferenceArray internal_face_edges = m->CreateTag("INTERNAL_FACE_EDGES",DATA_REFERENCE,NODE,NODE,4);
599
			//this marker helps detect edges of current cell only
Kirill Terekhov's avatar
Kirill Terekhov committed
600
			MarkerType mark_cell_edges = m->CreateMarker();
601
			//this marker helps detect nodes hanging on edges of unrefined cell
Kirill Terekhov's avatar
Kirill Terekhov committed
602
			MarkerType mark_hanging_nodes = m->CreateMarker();
603
			//9.split all cells of the current schedule
604
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
605
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
606
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
607
				Cell c = m->CellByLocalID(it);
608
609
610
611
612
613
614
615
616
617
618
				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
619
					Node n = m->CreateNode(xyz);
620
					//set increased level for the new node
621
					level[n] = level[c]+1;
622
623
624
625
626
					//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
627
628
					ElementArray<Node> edge_nodes(m,2);
					ElementArray<Edge> edges_to_faces(m,cell_hanging_nodes.size());
629
630
631
632
633
634
					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
635
						edges_to_faces[kt] = m->CreateEdge(edge_nodes).first;
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
669
670
						//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
671
					ElementArray<Face> internal_faces(m,edge_hanging_nodes.size());
672
673
674
675
676
677
678
679
680
681
					//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
682
						internal_faces[kt] = m->CreateFace(ElementArray<Edge>(m,face_edges.begin(),face_edges.end())).first;
683
684
685
686
687
						//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
688
689
690
691
692
693
694
					//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;
					//	
					//}
695
696
					//split the cell
					ElementArray<Cell> new_cells = Cell::SplitCell(c,internal_faces,0);
697
					std::sort(new_cells.begin(),new_cells.end(),Mesh::CentroidComparator(m));
698
					//retrive parent set
Kirill Terekhov's avatar
Kirill Terekhov committed
699
					ElementSet parent(m,parent_set[c]);
700
					//create set corresponding to old coarse cell
Kirill Terekhov's avatar
Kirill Terekhov committed
701
702
					Storage::real cnt[3];
					c.Centroid(cnt);
703
					std::stringstream set_name;
Kirill Terekhov's avatar
Kirill Terekhov committed
704
					//set_name << parent.GetName() << "_C" << c.GlobalID(); //rand may be unsafe
Kirill Terekhov's avatar
Kirill Terekhov committed
705
706
707
708
					if( parent == root )
						set_name << "r" << set_id[c];
					else
						set_name << parent.GetName() << "c" << set_id[c];
709
					//set_name << base64_encode_((unsigned char *)cnt,3*sizeof(double)/sizeof(unsigned char));
Kirill Terekhov's avatar
Kirill Terekhov committed
710
					
Kirill Terekhov's avatar
Kirill Terekhov committed
711
712
713
714
715
716
717
718
719
720
					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
721
						exit(-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
722
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
723
					
Kirill Terekhov's avatar
Kirill Terekhov committed
724
					ElementSet cell_set = m->CreateSetUnique(set_name.str()).first;
Kirill Terekhov's avatar
Kirill Terekhov committed
725
					cell_set->SetExchange(ElementSet::SYNC_ELEMENTS_ALL);
726
727
728
729
					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
730
						set_id[new_cells[kt]] = kt;
731
732
733
734
						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
735
736
737
738
739
740
741
742
743
744
					/*
					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() )
745
					parent.AddChild(cell_set);
Kirill Terekhov's avatar
Kirill Terekhov committed
746
					//else assert(cell_set->GetParent() == parent);
747
748
749
750
					//increment number of refined cells
					ret++;
				}
			}
751
			EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
752
753
754
			m->ReleaseMarker(mark_hanging_nodes);
			m->ReleaseMarker(mark_cell_edges);
			m->DeleteTag(internal_face_edges);
755
756
757
			//10.jump to later schedule, and go to 7.
			schedule_counter--;
		}
SilverLife's avatar
SilverLife committed
758

759
		//free created tag
Kirill Terekhov's avatar
Kirill Terekhov committed
760
		m->DeleteTag(indicator,FACE|EDGE);
SilverLife's avatar
SilverLife committed
761
762


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

Kirill Terekhov's avatar
Kirill Terekhov committed
787
788
789
790
791
792
        //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
793
    	//m->ResolveSets();
Kirill Terekhov's avatar
Kirill Terekhov committed
794
		
Kirill Terekhov's avatar
Kirill Terekhov committed
795
		//restore face orientation
Kirill Terekhov's avatar
Kirill Terekhov committed
796
		//BUG: bad orientation not fixed automatically
797
798
		
		/*
799
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
800
801
802
803
804
805
806
807
		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;
808
809
		if( nfixed ) REPORT_STR(rank << " fixed " << nfixed << " faces");
		EXIT_BLOCK();
810
811
		 */
		
Kirill Terekhov's avatar
Kirill Terekhov committed
812
813
814
815
816
817
818
819
820
		/*
		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
821
822
		//std::cout << rank << " check parent_set at end" << std::endl;
		//CheckParentSet();
823
824

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

    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;
    }
847
	
848
849
    void AdaptiveMesh::SynchronizeIndicated(TagInteger& indicator)
    {
850
        if (m->GetProcessorsNumber() == 1) return;
851
        ENTER_FUNC();
852
        int rank = m->GetProcessorRank();
853
854

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

                p++;
            }
        }
874
		EXIT_BLOCK();
875
876
        m->ReduceData(tag_indicated,ESET,0,OperationMin);
        m->ExchangeData(tag_indicated,ESET,0);
877
		ENTER_BLOCK();
878
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
879
880
881
882
883
884
885
886
887
            if (it->Integer(tag_indicated) == 0)
            {
                ElementSet::iterator p = it->Begin();
                while(p != it->End())
                {
                    p->Integer(indicator) = 0;
                    p++;
                }
            }
888
		EXIT_BLOCK();
889
        EXIT_FUNC();
890
891
    }

892
893
	bool AdaptiveMesh::Coarse(TagInteger & indicator)
	{
894
		std::string file;
895
896
		static int fi = 0;
        ENTER_FUNC();
Kirill Terekhov's avatar
Kirill Terekhov committed
897
        //return false;
898
899
900
901
		static int call_counter = 0;
		//return number of coarsened cells
		int ret = 0;
		//initialize tree structure
902
		ENTER_BLOCK();
903
		PrepareSet();
904
905
		EXIT_BLOCK();

906
907
		//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
908
		
909
		ENTER_BLOCK();
910
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
911
		m->ExchangeData(parent_set,CELL,0);
912
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
913
		m->ResolveSets();
914
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
915
		m->ExchangeData(parent_set,CELL,0);
916
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
917
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
918
919
		//m->CheckCentroids(__FILE__,__LINE__);
		//CheckParentSet();
920
921
		EXIT_BLOCK();

922
923
924
925
926
927
928
		//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
929
		
930
		
931
932
		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
933
934

		ENTER_BLOCK();
935
		//TagInteger coarsened = CreateTag("COARSENED",DATA_INTEGER,CELL,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
936
		TagInteger coarse_indicator = m->CreateTag("COARSE_INDICATOR",DATA_INTEGER,EDGE,NONE,1); //used to find on fine cells indicator on coarse cells
937
		//0. Extend indicator for sets, edges and faces
Kirill Terekhov's avatar
Kirill Terekhov committed
938
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
939
940
941
942
		//int mi = 0;
		//file = "w"+std::to_string(mi++)+".pvtk";
		//m->Save(file);
		//std::cout << "Save " << file << " " << __FILE__ << ":" << __LINE__ << std::endl;
943
944
		while(scheduled || unscheduled)
		{
945
946
			REPORT_VAL("scheduled",scheduled);
			REPORT_VAL("unscheduled",unscheduled);
947
948
949
950
951
952
			// 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
953
			m->ExchangeData(indicator,CELL,0);
954
955
956
			//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
957
			ENTER_BLOCK();
958
959
960
961
			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
962
				for(Storage::integer it = 0; it < m->LastLocalID(etype); ++it) if( m->isValidElement(etype,it) )
963
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
964
					Element e = m->ElementByLocalID(etype,it);
965
966
967
968
					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]]);
969
					//if (indicator[e] == INT_MAX) cout << rank << ": " << ElementTypeName(e.GetElementType()) << e.GlobalID() << endl;