amesh.cpp 49.7 KB
Newer Older
1
#include "amesh.h"
SilverLife's avatar
SilverLife committed
2
#include <iomanip>
3
#include <set>
Kirill Terekhov's avatar
Kirill Terekhov committed
4
#include "../../Source/Misc/base64.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
5
//using namespace std;
6

Kirill Terekhov's avatar
Kirill Terekhov committed
7
8
//from inmost
//std::string base64_encode(unsigned char const* buf, unsigned int bufLen);
Kirill Terekhov's avatar
Kirill Terekhov committed
9
/*
Kirill Terekhov's avatar
Kirill Terekhov committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
	std::string base64_encode_(unsigned char const* buf, unsigned int bufLen)
	{
		static const std::string base64_chars =
		"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
		"abcdefghijklmnopqrstuvwxyz"
		"0123456789+/";
		std::string ret;
		int i = 0;
		int j = 0;
		unsigned char char_array_3[3];
		unsigned char char_array_4[4];
		
		while (bufLen--)
		{
			char_array_3[i++] = *(buf++);
			if (i == 3)
			{
				char_array_4[0] = (char_array_3[0] & 0xfc) >> 2;
				char_array_4[1] = ((char_array_3[0] & 0x03) << 4) + ((char_array_3[1] & 0xf0) >> 4);
				char_array_4[2] = ((char_array_3[1] & 0x0f) << 2) + ((char_array_3[2] & 0xc0) >> 6);
				char_array_4[3] = char_array_3[2] & 0x3f;
				for(i = 0; (i <4) ; i++)
					ret += base64_chars[char_array_4[i]];
				i = 0;
			}
		}
		
		if (i)
		{
			for(j = i; j < 3; j++)
				char_array_3[j] = '\0';
	  
			char_array_4[0] = (char_array_3[0] & 0xfc) >> 2;
			char_array_4[1] = ((char_array_3[0] & 0x03) << 4) + ((char_array_3[1] & 0xf0) >> 4);
			char_array_4[2] = ((char_array_3[1] & 0x0f) << 2) + ((char_array_3[2] & 0xc0) >> 6);
			char_array_4[3] = char_array_3[2] & 0x3f;
	  
			for (j = 0; (j < i + 1); j++)
				ret += base64_chars[char_array_4[j]];
	  
			while((i++ < 3))
				ret += '=';
		}
		
		return ret;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
56
 */
57
58
59
/// todo:
/// 1. coarsment
/// 2. strategy for faces/edges with faults
Kirill Terekhov's avatar
Kirill Terekhov committed
60
/// 3. geom model supportbn
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
87
/// 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
88
89


Kirill Terekhov's avatar
Kirill Terekhov committed
90
    void AdaptiveMesh::PrintSetLocal(std::string offset, ElementSet it, std::stringstream& ss)
SilverLife's avatar
SilverLife committed
91
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
92
93
        std::stringstream ss1;
        ss1 << offset << rank << ": Set : " << std::setw(5) << it.GetName() << " ";
SilverLife's avatar
SilverLife committed
94
95
        for (int i = ss1.str().length(); i < 23; i++) ss1 << " ";
        ss << ss1.str();
Kirill Terekhov's avatar
Kirill Terekhov committed
96
        ss << std::setw(6);
SilverLife's avatar
SilverLife committed
97
98
99
100
101
        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";

102
        ss << " tag_owner (" << it.IntegerDF(m->OwnerTag()) << ")";
SilverLife's avatar
SilverLife committed
103
104
105

        //ss << "   level (" << level[it.self()] << ")  ";
        ss << " tag_processors (";
Kirill Terekhov's avatar
Kirill Terekhov committed
106
        std::stringstream ss2;
107
        Storage::integer_array arr = it.IntegerArrayDV(m->ProcessorsTag());
SilverLife's avatar
SilverLife committed
108
109
        for (int i = 0; i < arr.size(); i++)
            ss2 << arr[i] << " ";
Kirill Terekhov's avatar
Kirill Terekhov committed
110
        ss << std::setw(5) << ss2.str() <<")";
SilverLife's avatar
SilverLife committed
111
112
113
114
115
116
117
118
        

        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
119
            std::string type = "unknw";
SilverLife's avatar
SilverLife committed
120
121
122
123
            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
124
            ss << type << "-" << std::setw(2) << p->GlobalID() << " ";
SilverLife's avatar
SilverLife committed
125
126
            p++;
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
127
        ss << std::endl;
SilverLife's avatar
SilverLife committed
128
129
130
131
132
133
134
135
136
137

        for(ElementSet child = it.GetChild(); child.isValid(); child = child.GetSibling())
        {
            PrintSetLocal(offset + "   ",child,ss);
        }
    }


    void AdaptiveMesh::PrintSet()
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
138
        std::stringstream ss;
139
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
SilverLife's avatar
SilverLife committed
140
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
141
            //if (it->HaveParent()) continue;
142
            PrintSetLocal("",ElementSet(m,*it),ss);
SilverLife's avatar
SilverLife committed
143
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
144
        std::cout << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
145
146
    }

Kirill Terekhov's avatar
Kirill Terekhov committed
147
    void PrintRefs(std::ostream& os, Storage::reference_array refs)
148
149
150
    {
        for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
151
            std::string type = "unknw";
152
153
154
155
156
157
158
            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
159
            os << "(" << xyz[0] << "," << xyz[1] << "," << xyz[2] <<")" <<  std::endl;
160
161
162
163
164
        }


    }

Kirill Terekhov's avatar
Kirill Terekhov committed
165
    void AdaptiveMesh::PrintMesh(std::ostream& os, int cell, int face, int edge, int node)
SilverLife's avatar
SilverLife committed
166
167
    {
        if (cell + face + edge + node == 0) return;
Kirill Terekhov's avatar
Kirill Terekhov committed
168
169
        std::stringstream ss;
        ss << "================= " << rank << " =====================" << std::endl;
SilverLife's avatar
SilverLife committed
170
171
        if (cell)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
172
            ss << "Cells: " << m->NumberOfCells() <<  std::endl;
173
            for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) 
SilverLife's avatar
SilverLife committed
174
            {
175
                ss << rank << ": " << it->GlobalID() << " - " << it->LocalID() << " - ";            
SilverLife's avatar
SilverLife committed
176
177
178
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";
179
180
181
182
183
184
185
186
187
188
                
                /*
                ss << " (";
                Storage::integer_array arr = it->IntegerArrayDV(tag_processors);
                for (int i = 0; i < arr.size(); i++) ss << arr[i] << " ";
                ss << ") ";
*/
                /*
                Storage::reference_array refs = hanging_nodes[it->self()];
                if (refs.size() > 0)
SilverLife's avatar
SilverLife committed
189
                {
Kirill Terekhov's avatar
Kirill Terekhov committed
190
                    ss << std::endl << "   Hanging nodes: ";
191
                    PrintRefs(ss,refs);
SilverLife's avatar
SilverLife committed
192
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
193
                ss << std::endl;
194
195
196
197
198
199
                ss << "   ParentSet: " << ElementSet(this,parent_set[*it]).GetName();
                */
/*
                orage::reference_array refs = ref_tag[it->self()];
                PrintRefs(refs);
                */
Kirill Terekhov's avatar
Kirill Terekhov committed
200
                    ss << std::endl;
SilverLife's avatar
SilverLife committed
201
202
203
204
205
            }
        }

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

215
216
217

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

220
                ss << "  " << m->GetMarker(*it,m->NewMarker());
221
222
223
224

                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
225
                    ss << std::setw(2) << node->GlobalID() << " ";
226
227
228
                ss << ")";

                /*
SilverLife's avatar
SilverLife committed
229
230
231
232
233
234
235
236
237
238
239
240
                Storage::reference_array refs = ref_tag[it->self()];
                if (refs.size() > 0) ss << ". Ref: ";
                for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
                {
                    string type = "unknw";
                    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";
                    ss << "(" << type << "," << refs[i]->GlobalID() << ") ";
                }

Kirill Terekhov's avatar
Kirill Terekhov committed
241
                ss << std::endl;
242
243
244
245
246

                {
                    Storage::reference_array refs = hanging_nodes[it->self()];
                    if (refs.size() > 0)
                    {
Kirill Terekhov's avatar
Kirill Terekhov committed
247
                        ss << std::endl << "   Hanging nodes: ";
248
249
                        PrintRefs(ss,refs);
                    }
Kirill Terekhov's avatar
Kirill Terekhov committed
250
                    ss << std::endl;
251
                }
252
                */
Kirill Terekhov's avatar
Kirill Terekhov committed
253
                ss << std::endl;
SilverLife's avatar
SilverLife committed
254
255
256
257
258
            }
        }

        if (edge)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
259
            ss << "Edges: " << m->NumberOfEdges() <<  std::endl;
260
            for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it) 
SilverLife's avatar
SilverLife committed
261
262
263
264
265
            {
                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
266
				/*
SilverLife's avatar
SilverLife committed
267
268
269
270
                Storage::reference_array refs = ref_tag[it->self()];
                if (refs.size() > 0) ss << ". Ref: ";
                for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
                {
Kirill Terekhov's avatar
Kirill Terekhov committed
271
                    std::string type = "unknw";
SilverLife's avatar
SilverLife committed
272
273
274
275
276
277
                    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";
                    ss << "(" << type << "," << refs[i]->GlobalID() << ") ";
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
278
279
                */
                ss << std::endl;
SilverLife's avatar
SilverLife committed
280
281
282
283
284
            }
        }

        if (node)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
285
            ss << "Nodes:" << std::endl;
286
            for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) 
SilverLife's avatar
SilverLife committed
287
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
288
289
                ss << rank << ": " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
290
291
292
293
294
                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
295
					/*
SilverLife's avatar
SilverLife committed
296
297
298
299
                    Storage::reference_array refs = ref_tag[it->self()];
                    if (refs.size() > 0) ss << ". Ref: ";
                    for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
                    {
Kirill Terekhov's avatar
Kirill Terekhov committed
300
                        std::string type = "unknw";
SilverLife's avatar
SilverLife committed
301
                        if (refs[i].GetElementType() == CELL) type = "cell";
302
303
						if (refs[i].GetElementType() == FACE) type = "face";
						if (refs[i].GetElementType() == EDGE) type = "edge";
SilverLife's avatar
SilverLife committed
304
305
306
307
                        if (refs[i].GetElementType() == NODE) type = "node";
                        ss << "(" << type << "," << refs[i]->GlobalID() << ") ";
                        
                    }
Kirill Terekhov's avatar
Kirill Terekhov committed
308
					*/
309
                    ss << "  " << m->GetMarker(*it,m->NewMarker());
310
311

                    ss << "(" << 
Kirill Terekhov's avatar
Kirill Terekhov committed
312
313
314
                        std::setw(3) << it->RealArray(m->CoordsTag())[0] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[1] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[2] << ")";
315

SilverLife's avatar
SilverLife committed
316
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
317
                ss << std::endl;
SilverLife's avatar
SilverLife committed
318
319
320
            }
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
321
322
        ss << "=========================================" << std::endl;
        os << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
323
324
325
326
327
328
329
    }

    void AdaptiveMesh::UpdateStatus()
    {

        for(ElementType mask = CELL; mask >= NODE; mask = PrevElementType(mask))
        {
330
            for(Mesh::iteratorElement it = m->BeginElement(mask); it != m->EndElement(); it++)
SilverLife's avatar
SilverLife committed
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
            {
                int stat = 0;
                if (it->GetStatus() == Element::Shared) stat = 1;
                else if (it->GetStatus() == Element::Ghost)  stat = 2;

                tag_status[it->self()] = stat;
            }
        }
    }
    
    void AdaptiveMesh::PrintSet(ElementSet set, std::string offset)
    {
        std::cout << offset << "Set: " << std::endl;

		ElementSet::iterator it = set.Begin();
		while(it != set.End())
        {
348
		    ElementSet parent(m,parent_set[*it]);
SilverLife's avatar
SilverLife committed
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
            std::cout << offset << it->GlobalID() << " - " << level[*it] << " : ";
		    
            ElementSet::iterator p = parent.Begin();
    		while(p != parent.End())
            {
                std::cout << p->GlobalID() << " ";
                p++;
            }
            std::cout << std::endl;
            it++;
        }

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

    }

368
369
370
    void AdaptiveMesh::SynchronizeSet(ElementSet set)
    {
    #ifdef USE_MPI
371
372
        int size = m->GetProcessorsNumber();
        int rank = m->GetProcessorRank();
373
374
        for (int i = 0; i < size; i++)
        {
375
376
            set.IntegerArray(m->SendtoTag()).push_back(i);
            m->ExchangeMarked();
377
378
379
380
        }
    #endif
    }

SilverLife's avatar
SilverLife committed
381
382
    void AdaptiveMesh::Test()
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
383
        std::cout << rank << ": ================" << std::endl;
SilverLife's avatar
SilverLife committed
384
385
        PrintSet(root,"");
    }
386
387
388
389
	
	
	void AdaptiveMesh::ClearData()
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
390
		set_id        = m->DeleteTag(set_id);
Kirill Terekhov's avatar
Kirill Terekhov committed
391
392
393
		level         = m->DeleteTag(level);
		hanging_nodes = m->DeleteTag(hanging_nodes);
		parent_set    = m->DeleteTag(parent_set);
394
395
396
397
398
399
400
401
		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
402
			root = m->GetSet("ROOT_SET");
403
404
			if( root == InvalidElement() )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
405
				root = m->CreateSetUnique("ROOT_SET").first;
Kirill Terekhov's avatar
Kirill Terekhov committed
406
				root.SetExchange(ElementSet::SYNC_ELEMENTS_SHARED);
407
				level[root] = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
408
				for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
409
410
411
412
				{
					root.PutElement(it->self());
					parent_set[it->self()] = root.GetHandle();
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
413
				m->Enumerate(CELL,set_id);
414
			}
SilverLife's avatar
SilverLife committed
415
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
416
		if( !m->HaveGlobalID(CELL) ) m->AssignGlobalID(CELL); //for unique set names
Kirill Terekhov's avatar
Kirill Terekhov committed
417
		m->ResolveSets();
418
419
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
420
	AdaptiveMesh::AdaptiveMesh(Mesh & _m) : m(&_m)
421
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
422
		model = NULL;
423
		//create a tag that stores maximal refinement level of each element
Kirill Terekhov's avatar
Kirill Terekhov committed
424
		level = m->CreateTag("REFINEMENT_LEVEL",DATA_INTEGER,CELL|FACE|EDGE|NODE|ESET,NONE,1);
425
		tag_status = m->CreateTag("TAG_STATUS",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
426
		set_id = m->CreateTag("SET_ID",DATA_INTEGER,CELL,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
427
428
		//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);
429
		//create a tag that stores links to all the hanging nodes of the cell
Kirill Terekhov's avatar
Kirill Terekhov committed
430
		hanging_nodes = m->CreateTag("HANGING_NODES",DATA_REFERENCE,CELL|FACE,NONE);
431
		//create a tag that stores links to sets
Kirill Terekhov's avatar
Kirill Terekhov committed
432
		parent_set = m->CreateTag("PARENT_SET",DATA_REFERENCE,CELL,NONE,1);
433
434
	    size = m->GetProcessorsNumber();
    	rank = m->GetProcessorRank();
435
436
437
438
439
440
441
442
	}
	
	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
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
	void AdaptiveMesh::CheckParentSet()
	{
		std::cout << rank << " enter " << __FUNCTION__ << std::endl;
		int err = 0;
		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()) << " " << level[*it] << std::endl;
				err++;
			}
			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()) << " " << level[*it] << std::endl;
				err++;
			}
		}
		err = m->Integrate(err);
		if( err ) 
		{
			std::cout << rank << " error in " << __FUNCTION__ << std::endl;
			exit(-1);
		}
		std::cout << rank << " exit " << __FUNCTION__ << std::endl;
		
	}
	
470
471
	bool AdaptiveMesh::Refine(TagInteger & indicator)
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
472
        std::cout << rank << " enter " << __FUNCTION__ << std::endl;
473
474
475
		static int call_counter = 0;
		int ret = 0; //return number of refined cells
		//initialize tree structure
Kirill Terekhov's avatar
Kirill Terekhov committed
476
		m->CheckCentroids(__FILE__,__LINE__);
477
		PrepareSet();
Kirill Terekhov's avatar
Kirill Terekhov committed
478
479
480
481
482
483
484
485
486
487
488
		
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(parent_set,CELL,0);
		m->CheckCentroids(__FILE__,__LINE__);
		m->ResolveSets();
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(parent_set,CELL,0);
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
		m->CheckCentroids(__FILE__,__LINE__);
		
Kirill Terekhov's avatar
Kirill Terekhov committed
489
		CheckParentSet();
490
491
492
		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
		//0. Extend indicator for edges and faces
Kirill Terekhov's avatar
Kirill Terekhov committed
493
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
494
495
496
		while(scheduled)
		{
			//1.Communicate indicator - it may be not synced
Kirill Terekhov's avatar
Kirill Terekhov committed
497
			m->ExchangeData(indicator,CELL,0);
498
499
			//2.Propogate indicator down to the faces,edges
			//  select schedule for them
Kirill Terekhov's avatar
Kirill Terekhov committed
500
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
501
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
502
				Cell c = m->CellByLocalID(it);
503
504
505
506
507
508
509
510
511
512
513
				if( indicator[c] == schedule_counter )
				{
					ElementArray<Element> adj = c.getAdjElements(FACE|EDGE);
					for(ElementArray<Element>::size_type kt = 0; kt < adj.size(); ++kt)
					{
						if( level[adj[kt]] == level[c] ) //do not schedule finer or coarser elements
							indicator[adj[kt]] = schedule_counter; //refine together with current cell
					}
				}
			}
			//3.Communicate indicator on faces and edges
Kirill Terekhov's avatar
Kirill Terekhov committed
514
			m->ExchangeData(indicator,FACE|EDGE,0);
515
516
517
518
			//4.Check for each cell if there is
			//  any hanging node with adjacent in a need to refine,
			//  schedule for refinement earlier.
			scheduled = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
519
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
520
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
521
				Cell c = m->CellByLocalID(it);
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
				//already scheduled cells may be required to be refined first
				//if( indicator[c] == 0 ) //some optimization
				{
					bool scheduled_c = false;
					//any finer level edge is scheduled to be refined first
					ElementArray<Edge> edges = c->getEdges();
					for(ElementArray<Edge>::size_type kt = 0; kt < edges.size() && !scheduled_c; ++kt)
					{
						//if a finer edge is scheduled
						//then this cell should be refined first
						if( indicator[edges[kt]] != 0 &&
							level[edges[kt]] > level[c] &&
							indicator[edges[kt]] >= indicator[c] )
						{
							indicator[c] = schedule_counter+1;
							scheduled++;
							scheduled_c = true;
						}
					}
				}
			}
			//5.Go back to 1 until no new elements scheduled
Kirill Terekhov's avatar
Kirill Terekhov committed
544
			scheduled = m->Integrate(scheduled);
545
546
			if( scheduled ) schedule_counter++;
		}
547
		m->ExchangeData(indicator,CELL | FACE | EDGE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
548
		//m->Save("indicator.pmf");
549
		//6.Refine
Kirill Terekhov's avatar
Kirill Terekhov committed
550
		m->BeginModification();
551
552
553
554
		while(schedule_counter)
		{
			Storage::real xyz[3] = {0,0,0};
			//7.split all edges of the current schedule
Kirill Terekhov's avatar
Kirill Terekhov committed
555
			for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
556
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
557
				Edge e = m->EdgeByLocalID(it);
558
559
560
561
562
				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
563
					for(Storage::integer d = 0; d < m->GetDimensions(); ++d)
564
565
566
						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
567
					Node n = m->CreateNode(xyz);
568
					//set increased level for new node
569
					level[n] = level[e.getBeg()] = level[e.getEnd()] = level[e]+1;
570
571
572
573
					//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
574
					ElementArray<Edge> new_edges = Edge::SplitEdge(e,ElementArray<Node>(m,1,n.GetHandle()),0);
575
576
577
578
579
					//set increased level for new edges
					level[new_edges[0]] = level[new_edges[1]] = level[e]+1;
				}
			}
			//8.split all faces of the current schedule, using hanging nodes on edges
Kirill Terekhov's avatar
Kirill Terekhov committed
580
			for(Storage::integer it = 0; it < m->FaceLastLocalID(); ++it) if( m->isValidFace(it) )
581
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
582
				Face f = m->FaceByLocalID(it);
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
				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
598
					Node n = m->CreateNode(xyz);
599
					//set increased level for the new node
600
					level[n] = level[f]+1;
601
602
603
					//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
604
605
					ElementArray<Node> edge_nodes(m,2); //to create new edges
					ElementArray<Edge> hanging_edges(m,face_hanging_nodes.size());
606
607
608
609
					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
610
						hanging_edges[kt] = m->CreateEdge(edge_nodes).first;
611
612
613
614
615
616
617
618
619
620
621
						//set increased level for new edges
						level[hanging_edges[kt]] = level[f]+1;
					}
					//split the face by these edges
					ElementArray<Face> new_faces = Face::SplitFace(f,hanging_edges,0);
					//set increased level to new faces
					for(ElementArray<Face>::size_type kt = 0; kt < new_faces.size(); ++kt)
						level[new_faces[kt]] = level[f]+1;
				}
			}
			//this tag helps recreate internal face
Kirill Terekhov's avatar
Kirill Terekhov committed
622
			TagReferenceArray internal_face_edges = m->CreateTag("INTERNAL_FACE_EDGES",DATA_REFERENCE,NODE,NODE,4);
623
			//this marker helps detect edges of current cell only
Kirill Terekhov's avatar
Kirill Terekhov committed
624
			MarkerType mark_cell_edges = m->CreateMarker();
625
			//this marker helps detect nodes hanging on edges of unrefined cell
Kirill Terekhov's avatar
Kirill Terekhov committed
626
			MarkerType mark_hanging_nodes = m->CreateMarker();
627
			//9.split all cells of the current schedule
Kirill Terekhov's avatar
Kirill Terekhov committed
628
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
629
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
630
				Cell c = m->CellByLocalID(it);
631
632
633
634
635
636
637
638
639
640
641
				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
642
					Node n = m->CreateNode(xyz);
643
					//set increased level for the new node
644
					level[n] = level[c]+1;
645
646
647
648
649
					//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
650
651
					ElementArray<Node> edge_nodes(m,2);
					ElementArray<Edge> edges_to_faces(m,cell_hanging_nodes.size());
652
653
654
655
656
657
					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
658
						edges_to_faces[kt] = m->CreateEdge(edge_nodes).first;
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
						//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
694
					ElementArray<Face> internal_faces(m,edge_hanging_nodes.size());
695
696
697
698
699
700
701
702
703
704
					//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
705
						internal_faces[kt] = m->CreateFace(ElementArray<Edge>(m,face_edges.begin(),face_edges.end())).first;
706
707
708
709
710
						//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
711
712
713
714
715
716
717
					//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;
					//	
					//}
718
719
720
					//split the cell
					ElementArray<Cell> new_cells = Cell::SplitCell(c,internal_faces,0);
					//retrive parent set
Kirill Terekhov's avatar
Kirill Terekhov committed
721
					ElementSet parent(m,parent_set[c]);
722
					//create set corresponding to old coarse cell
Kirill Terekhov's avatar
Kirill Terekhov committed
723
724
					Storage::real cnt[3];
					c.Centroid(cnt);
725
					std::stringstream set_name;
Kirill Terekhov's avatar
Kirill Terekhov committed
726
					//set_name << parent.GetName() << "_C" << c.GlobalID(); //rand may be unsafe
Kirill Terekhov's avatar
Kirill Terekhov committed
727
728
729
730
					if( parent == root )
						set_name << "r" << set_id[c];
					else
						set_name << parent.GetName() << "c" << set_id[c];
Kirill Terekhov's avatar
Kirill Terekhov committed
731
					set_name << base64_encode_((unsigned char *)cnt,3*sizeof(double)/sizeof(unsigned char));
Kirill Terekhov's avatar
Kirill Terekhov committed
732
					
Kirill Terekhov's avatar
Kirill Terekhov committed
733
734
735
736
737
738
739
740
741
742
					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
743
						exit(-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
744
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
745
					
Kirill Terekhov's avatar
Kirill Terekhov committed
746
					ElementSet cell_set = m->CreateSetUnique(set_name.str()).first;
Kirill Terekhov's avatar
Kirill Terekhov committed
747
					cell_set->SetExchange(ElementSet::SYNC_ELEMENTS_ALL);
748
749
750
751
					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
752
						set_id[new_cells[kt]] = kt;
753
754
755
756
						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
757
758
759
760
761
762
763
764
765
766
					/*
					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() )
767
					parent.AddChild(cell_set);
Kirill Terekhov's avatar
Kirill Terekhov committed
768
					//else assert(cell_set->GetParent() == parent);
769
770
771
772
					//increment number of refined cells
					ret++;
				}
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
773
774
775
			m->ReleaseMarker(mark_hanging_nodes);
			m->ReleaseMarker(mark_cell_edges);
			m->DeleteTag(internal_face_edges);
776
777
778
			//10.jump to later schedule, and go to 7.
			schedule_counter--;
		}
SilverLife's avatar
SilverLife committed
779

780
		//free created tag
Kirill Terekhov's avatar
Kirill Terekhov committed
781
		m->DeleteTag(indicator,FACE|EDGE);
SilverLife's avatar
SilverLife committed
782

783
		/*
784
785
        MarkerType marker_new = m->CreateMarker();
        for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) 
SilverLife's avatar
SilverLife committed
786
        {
787
            if (it->GetMarker(m->NewMarker()) == false) continue;
788
789
790
791
            it->SetMarker(marker_new);

        }

792
        for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) 
793
        {
794
            if (it->GetMarker(m->NewMarker()) == false) continue;
795
796
            it->SetMarker(marker_new);
        }
SilverLife's avatar
SilverLife committed
797

798
        for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it) 
799
        {
800
            if (it->GetMarker(m->NewMarker()) == false) continue;
801
            it->SetMarker(marker_new);
SilverLife's avatar
SilverLife committed
802
        }
803

804
        for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) 
SilverLife's avatar
SilverLife committed
805
        {
806
            if (it->GetMarker(m->NewMarker()) == false) continue;
807
            it->SetMarker(marker_new);
SilverLife's avatar
SilverLife committed
808
        }
809
        */
SilverLife's avatar
SilverLife committed
810
811


812
		//11. Restore parallel connectivity, global ids
Kirill Terekhov's avatar
Kirill Terekhov committed
813
		/*
Kirill Terekhov's avatar
Kirill Terekhov committed
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
		{
			std::fstream fout;
			fout.open("out"+std::to_string(m->GetProcessorRank())+".txt",std::ios::out);
			for(ElementType etype = NODE; etype <= ESET; etype = NextElementType(etype) )
			{
				int cnt = 0;
				for(Mesh::iteratorElement it = m->BeginElement(etype); it != m->EndElement(); ++it)
				{
					if( it->New() )
					{
						fout << "new element " << ElementTypeName(etype) << ":" << it->LocalID() <<std::endl;//<< " gid " << it->GlobalID() << std::endl;
						cnt++;
					}
				}
				fout << ElementTypeName(etype) << " count: " << cnt << std::endl;
			}
			
			fout.close();
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
833
		 */
SilverLife's avatar
SilverLife committed
834
        //if (call_counter == 0)
Kirill Terekhov's avatar
Kirill Terekhov committed
835
		m->ResolveModification();
Kirill Terekhov's avatar
Kirill Terekhov committed
836
		//m->SynchronizeMarker(m->NewMarker(),CELL|FACE|EDGE|NODE,SYNC_BIT_OR);
837
        //ExchangeGhost(3,NODE); // Construct Ghost cells in 2 layers connected via nodes
838
839
		//12. Let the user update their data
		//todo: call back function for New() cells
Kirill Terekhov's avatar
Kirill Terekhov committed
840
		if( model ) model->Adaptation(*m);
841
		//13. Delete old elements of the mesh
Kirill Terekhov's avatar
Kirill Terekhov committed
842
		m->ApplyModification();
Kirill Terekhov's avatar
Kirill Terekhov committed
843
		
Kirill Terekhov's avatar
Kirill Terekhov committed
844
		//m->ExchangeGhost(1,NODE,m->NewMarker());
845
		//14. Done
Kirill Terekhov's avatar
Kirill Terekhov committed
846
        //cout << rank << ": Before end " << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
847
		m->EndModification();
Kirill Terekhov's avatar
Kirill Terekhov committed
848
		
Kirill Terekhov's avatar
Kirill Terekhov committed
849
850
851
852
		static int fi = 0;
		m->Save("refine"+std::to_string(fi)+".pvtk");
		std::cout << "Save refine"+std::to_string(fi)+".pvtk" << std::endl;
		fi++;
Kirill Terekhov's avatar
Kirill Terekhov committed
853
		
854
		//ExchangeData(hanging_nodes,CELL | FACE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
855
        //m->ResolveSets();
856

Kirill Terekhov's avatar
Kirill Terekhov committed
857
858
859
860
861
862
        //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
863
    	//m->ResolveSets();
Kirill Terekhov's avatar
Kirill Terekhov committed
864
		
Kirill Terekhov's avatar
Kirill Terekhov committed
865
		//restore face orientation
Kirill Terekhov's avatar
Kirill Terekhov committed
866
867
868
869
870
871
872
873
874
		//BUG: bad orientation not fixed automatically
		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;
Kirill Terekhov's avatar
Kirill Terekhov committed
875
876
877
878
879
880
881
882
883
884
		if( nfixed ) std::cout << rank << " fixed " << nfixed << " faces" << std::endl;
		/*
		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
885
886
		//std::cout << rank << " check parent_set at end" << std::endl;
		//CheckParentSet();
887
888

        //ExchangeData(hanging_nodes,CELL | FACE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
889
        //cout << rank << ": After end " << std::endl;
890
		//reorder element's data to free up space
Kirill Terekhov's avatar
Kirill Terekhov committed
891
		m->ReorderEmpty(CELL|FACE|EDGE|NODE);
892
893
		//return number of refined cells
		call_counter++;
Kirill Terekhov's avatar
Kirill Terekhov committed
894
        std::cout << rank << " exit " << __FUNCTION__ << std::endl;
895
896
		return ret != 0;
	}
897
898
899
900
901
902
903
904
905
906
907
908

    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;
    }
909
	
910
911
    void AdaptiveMesh::SynchronizeIndicated(TagInteger& indicator)
    {
912
        if (m->GetProcessorsNumber() == 1) return;
Kirill Terekhov's avatar
Kirill Terekhov committed
913
        std::cout << rank << " enter " << __FUNCTION__ << std::endl;
914
        int rank = m->GetProcessorRank();
915
916

        // Check all sets. All elements in sets must be indicated. At first we check indicator in local processor, and second integrate data
917
918
        TagInteger tag_indicated = m->CreateTag("INDICATED",DATA_INTEGER,ESET,NONE,1);
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
919
        {
920
            ElementSet set = ElementSet(m,*it);
921
922
923
924
925
926
927
            set.Integer(tag_indicated) = 1;
            ElementSet::iterator p = set.Begin();
            while(p != set.End())
            {
                if (indicator[*p] == 0)
                {
                    tag_indicated[set] = 0;
928
                    //std::cout << rank << ": Set " << set.GetName() << " not all indicated" << endl;
929
930
931
932
933
934
935
                    break;
                }

                p++;
            }
        }

936
937
        m->ReduceData(tag_indicated,ESET,0,OperationMin);
        m->ExchangeData(tag_indicated,ESET,0);
938

939
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
            if (it->Integer(tag_indicated) == 0)
            {
                ElementSet::iterator p = it->Begin();
                while(p != it->End())
                {
                    p->Integer(indicator) = 0;
                    p++;
                }
            }
        /*
                stringstream ss;
        for(Mesh::iteratorSet it = BeginSet(); it != EndSet(); ++it) 
        {
            ss << it->GetName() << " - " << it->Integer(tag_indicated) << endl;;
        }
955
        cout << rank << " Sets: \n" << ss.str() << endl;
956
        */
Kirill Terekhov's avatar
Kirill Terekhov committed
957
        std::cout << rank << " exit " << __FUNCTION__ << std::endl;
958
959
    }

960
961
	bool AdaptiveMesh::Coarse(TagInteger & indicator)
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
962
        std::cout << rank << " enter " << __FUNCTION__ << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
963
        //return false;
Kirill Terekhov's avatar
Kirill Terekhov committed
964
		m->CheckCentroids(__FILE__,__LINE__);
965
966
967
968
969
		static int call_counter = 0;
		//return number of coarsened cells
		int ret = 0;
		//initialize tree structure
		PrepareSet();
Kirill Terekhov's avatar
Kirill Terekhov committed
970
971
972
973
974
975
976
977
978
979
980
		
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(parent_set,CELL,0);
		m->CheckCentroids(__FILE__,__LINE__);
		m->ResolveSets();
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(parent_set,CELL,0);
		m->CheckCentroids(__FILE__,__LINE__);
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
		m->CheckCentroids(__FILE__,__LINE__);
		
Kirill Terekhov's avatar
Kirill Terekhov committed
981
982
		SynchronizeIndicated(indicator);
		CheckParentSet();
983
984
985
		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
		//TagInteger coarsened = CreateTag("COARSENED",DATA_INTEGER,CELL,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
986
		TagInteger coarse_indicator = m->CreateTag("COARSE_INDICATOR",DATA_INTEGER,EDGE,NONE,1); //used to find on fine cells indicator on coarse cells
987
		//0. Extend indicator for sets, edges and faces
Kirill Terekhov's avatar
Kirill Terekhov committed
988
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
989
990
991
992
993
994
995
996
		while(scheduled || unscheduled)
		{
			// rules
			// a) If there is adjacent finer edge that is not marked for coarsening
			// then this cell should not be coarsened
			// b) If there is adjacent coarser cell, then this cell should be coarsened
			// first
			//0.Communicate indicator - it may be not synced
Kirill Terekhov's avatar
Kirill Terekhov committed
997
			m->ExchangeData(indicator,CELL,0);
998
999
1000
1001
1002
1003
1004
			//1. Mark each adjacent face/edge for coarsement schedule
			// problem: should mark so that if every adjacent cell is coarsened
			// then adjacent face/edge are also coarsened
			for(ElementType etype = EDGE; etype <= FACE; etype = NextElementType(etype))
			{
				//for(Storage::integer it = 0; it < LastLocalID(etype); ++it) if( isValidElement(etype,it) )
				//	indicator[ElementByLocalID(etype,it)] = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1005
				for(Storage::integer it = 0; it < m->LastLocalID(etype); ++it) if( m->isValidElement(etype,it) )
1006
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
1007
					Element e = m->ElementByLocalID(etype,it);
1008
1009
1010
1011
					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]]);
1012
					//if (indicator[e] == INT_MAX) cout << rank << ": " << ElementTypeName(e.GetElementType()) << e.GlobalID() << endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
1013
					//assert(indicator[e] != INT_MAX);
1014
1015
1016
				}
			}
			//2.Communicate indicator on faces and edges
Kirill Terekhov's avatar
Kirill Terekhov committed
1017
			m->ReduceData(indicator,FACE|EDGE,0,ReduceMin);
1018
			m->ExchangeData(indicator,FACE|EDGE,0);
1019
1020
1021
			//3.If there is adjacent finer edge that are not marked for coarsening
			// then this cell should not be coarsened
			unscheduled = scheduled = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1022
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
1023
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1024
				Cell c = m->CellByLocalID(it);
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
				if( indicator[c] )
				{
					ElementArray<Edge> edges = c.getEdges();
					for(ElementArray<Edge>::size_type kt = 0; kt < edges.size(); ++kt)
					{
						if( level[edges[kt]] > level[c] && indicator[edges[kt]] == 0 )
						{
							indicator[c] = 0;
							unscheduled++;
						}
					}
				}
			}
			//4. Propogate coarsement info over set tree to detect valid coarsenings.
			// go down over sets, if set does not have children and all of the cells
			// of the set are marked for coarsening, then mark the set for coarsement
			// otherwise unmark.
			// Unmark all cells that are not to be coarsened
Kirill Terekhov's avatar
Kirill Terekhov committed
1043
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
1044
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1045
				Cell c = m->CellByLocalID(it);
1046
                //if (!isValidElement(c.GetHandle())) continue;
1047
1048
				if( indicator[c] )
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
1049
					ElementSet parent(m,parent_set[c]);
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
					//intermediate cell may not be coarsened
					//root set may not have coarsening cells
					if( parent.HaveChild() || !parent.HaveParent() )
					{
						indicator[c] = 0;
						unscheduled++;
					}
					else
					{
						Storage::integer schedule_first = 0;
						bool check = true;
						//check that all elements of the set are to be coarsened
						for(ElementSet::iterator it = parent.Begin(); it != parent.End(); ++it)
						{
							check &= (indicator[it->self()] != 0);
							schedule_first = std::max(schedule_first,indicator[it->self()]);
						}
						if(!check)
						{
							indicator[c] = 0;
							unscheduled++;
						}
						else if( indicator[c] != schedule_first )
						{
							indicator[c] = schedule_first;
							unscheduled++;
						}
					}
				}
			}
			//5.If there is an adjacent coarser element to be refined, then
			//   this one should be scheduled to be refined first
			//a) clean up coarse indicator tag
Kirill Terekhov's avatar
Kirill Terekhov committed
1083
1084
			for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
				coarse_indicator[m->EdgeByLocalID(it)] = 0;
1085
			//b) each cell mark it's finer edges with cell's schedule
Kirill Terekhov's avatar
Kirill Terekhov committed
1086
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
1087
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1088
				Cell c = m->CellByLocalID(it);
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
				if( indicator[c] )
				{
					ElementArray<Element> adj = c.getAdjElements(EDGE);
					for(ElementArray<Element>::size_type kt = 0; kt < adj.size(); ++kt)
					{
						if( level[adj[kt]] > level[c] ) //only finer edges
							coarse_indicator[adj[kt]] = std::max(coarse_indicator[adj[kt]],indicator[c]);
					}
				}
			}
			//c) data reduction to get maximum over mesh partition
Kirill Terekhov's avatar
Kirill Terekhov committed
1100
			m->ReduceData(coarse_indicator,EDGE,0,ReduceMax);
1101
			m->ExchangeData(coarse_indicator,EDGE,0);
1102
			//d) look from cells if any edge is coarsened earlier
Kirill Terekhov's avatar
Kirill Terekhov committed
1103
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
1104
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1105
				Cell c = m->CellByLocalID(it);
1106
1107
1108
1109
1110
1111
1112
1113
1114