amesh.cpp 64.9 KB
Newer Older
1
#include "amesh.h"
SilverLife's avatar
SilverLife committed
2
#include <iomanip>
3
#include <set>
4
5
6
7
8
9
10


#if defined(USE_PARALLEL_WRITE_TIME)
#define REPORT_MPI(x) {m->WriteTab(m->out_time) << "<MPI><![CDATA[" << #x << "]]></MPI>" << std::endl; x;}
#define REPORT_STR(x) {m->WriteTab(m->out_time) << "<TEXT><![CDATA[" << x << "]]></TEXT>" << std::endl;}
#define REPORT_VAL(str,x) {m->WriteTab(m->out_time) << "<VALUE name=\"" << str << "\"> <CONTENT><![CDATA[" << x << "]]></CONTENT> <CODE><![CDATA[" << #x << "]]></CODE></VALUE>" << std::endl;}
#define ENTER_FUNC() double all_time = Timer(); {m->WriteTab(m->out_time) << "<FUNCTION name=\"" << __FUNCTION__ << "\" id=\"func" << m->GetFuncID()++ << "\">" << std::endl; m->Enter();}
11
#define ENTER_BLOCK() { double btime = Timer(); m->WriteTab(m->out_time) << "<FUNCTION name=\"" << __FUNCTION__ << ":" << NameSlash(__FILE__) << ":" << __LINE__ << "\" id=\"func" << m->GetFuncID()++ << "\">" << std::endl; m->Enter();
12
13
14
15
16
17
18
#define EXIT_BLOCK() m->WriteTab(m->out_time) << "<TIME>" << Timer() - btime << "</TIME>" << std::endl; m->Exit(); m->WriteTab(m->out_time) << "</FUNCTION>" << std::endl;}
#define EXIT_FUNC() {m->WriteTab(m->out_time) << "<TIME>" << Timer() - all_time << "</TIME>" << std::endl; m->Exit(); m->WriteTab(m->out_time) << "</FUNCTION>" << std::endl;}
#define EXIT_FUNC_DIE() {m->WriteTab(m->out_time) << "<TIME>" << -1 << "</TIME>" << std::endl; m->Exit(); m->WriteTab(m->out_time) << "</FUNCTION>" << std::endl;}
#else
#define REPORT_MPI(x) x
#define REPORT_STR(x) {}
#define REPORT_VAL(str,x) {}
19
20
#define ENTER_BLOCK()
#define EXIT_BLOCK()
21
22
23
24
25
#define ENTER_FUNC() {}
#define EXIT_FUNC() {}
#define EXIT_FUNC_DIE()  {}
#endif

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

34

35
36
37
/// todo:
/// 1. coarsment
/// 2. strategy for faces/edges with faults
Kirill Terekhov's avatar
Kirill Terekhov committed
38
/// 3. geom model supportbn
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
/// 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));
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
66
67
68
69
70
71
72
73
74
75
	
	void ReduceUnion(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void) size;
		const INMOST_DATA_INTEGER_TYPE * idata = (const INMOST_DATA_INTEGER_TYPE *)data;
		Storage::integer_array odata = element->IntegerArray(tag);
		std::vector<int> tmp(size+odata.size());
		tmp.resize(std::set_union(odata.begin(),odata.end(),idata,idata+size,tmp.begin())-tmp.begin());
		odata.replace(odata.begin(),odata.end(),tmp.begin(),tmp.end());
	}
SilverLife's avatar
SilverLife committed
76

77
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
78
    void AdaptiveMesh::PrintSetLocal(std::string offset, ElementSet it, std::stringstream& ss)
SilverLife's avatar
SilverLife committed
79
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
80
81
        std::stringstream ss1;
        ss1 << offset << rank << ": Set : " << std::setw(5) << it.GetName() << " ";
SilverLife's avatar
SilverLife committed
82
83
        for (int i = ss1.str().length(); i < 23; i++) ss1 << " ";
        ss << ss1.str();
Kirill Terekhov's avatar
Kirill Terekhov committed
84
        ss << std::setw(6);
SilverLife's avatar
SilverLife committed
85
86
87
88
89
        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";

90
        ss << " tag_owner (" << it.IntegerDF(m->OwnerTag()) << ")";
SilverLife's avatar
SilverLife committed
91
92
93

        //ss << "   level (" << level[it.self()] << ")  ";
        ss << " tag_processors (";
Kirill Terekhov's avatar
Kirill Terekhov committed
94
        std::stringstream ss2;
95
        Storage::integer_array arr = it.IntegerArrayDV(m->ProcessorsTag());
SilverLife's avatar
SilverLife committed
96
97
        for (int i = 0; i < arr.size(); i++)
            ss2 << arr[i] << " ";
Kirill Terekhov's avatar
Kirill Terekhov committed
98
        ss << std::setw(5) << ss2.str() <<")";
SilverLife's avatar
SilverLife committed
99
100
101
102
103
104
105
106
        

        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
107
            std::string type = "unknw";
SilverLife's avatar
SilverLife committed
108
109
110
111
            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
112
            ss << type << "-" << std::setw(2) << p->GlobalID() << " ";
SilverLife's avatar
SilverLife committed
113
114
            p++;
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
115
        ss << std::endl;
SilverLife's avatar
SilverLife committed
116
117
118
119
120
121

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

124
	/*
SilverLife's avatar
SilverLife committed
125
126
    void AdaptiveMesh::PrintSet()
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
127
        std::stringstream ss;
128
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
SilverLife's avatar
SilverLife committed
129
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
130
            //if (it->HaveParent()) continue;
131
            PrintSetLocal("",ElementSet(m,*it),ss);
SilverLife's avatar
SilverLife committed
132
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
133
        std::cout << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
134
    }
135
	*/
SilverLife's avatar
SilverLife committed
136

137
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
138
    void PrintRefs(std::ostream& os, Storage::reference_array refs)
139
140
141
    {
        for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
142
            std::string type = "unknw";
143
144
145
146
147
148
149
            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
150
            os << "(" << xyz[0] << "," << xyz[1] << "," << xyz[2] <<")" <<  std::endl;
151
152
        }
    }
153
	*/
154

155
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
156
    void AdaptiveMesh::PrintMesh(std::ostream& os, int cell, int face, int edge, int node)
SilverLife's avatar
SilverLife committed
157
158
    {
        if (cell + face + edge + node == 0) return;
Kirill Terekhov's avatar
Kirill Terekhov committed
159
160
        std::stringstream ss;
        ss << "================= " << rank << " =====================" << std::endl;
SilverLife's avatar
SilverLife committed
161
162
        if (cell)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
163
            ss << "Cells: " << m->NumberOfCells() <<  std::endl;
164
            for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) 
SilverLife's avatar
SilverLife committed
165
            {
166
                ss << rank << ": " << it->GlobalID() << " - " << it->LocalID() << " - ";            
SilverLife's avatar
SilverLife committed
167
168
169
                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
170
                ss << std::endl;
SilverLife's avatar
SilverLife committed
171
172
173
174
175
            }
        }

        if (face)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
176
            ss << "Faces: " << m->NumberOfFaces() <<  std::endl;
177
            for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) 
SilverLife's avatar
SilverLife committed
178
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
179
180
                ss << rank << ": " << std::setw(2) << it->LocalID() << " " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
181
182
183
184
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";

185
186
187

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

190
                ss << "  " << m->GetMarker(*it,m->NewMarker());
191
192
193
194

                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
195
                    ss << std::setw(2) << node->GlobalID() << " ";
196
                ss << ")";
Kirill Terekhov's avatar
Kirill Terekhov committed
197
                ss << std::endl;
SilverLife's avatar
SilverLife committed
198
199
200
201
202
            }
        }

        if (edge)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
203
            ss << "Edges: " << m->NumberOfEdges() <<  std::endl;
204
            for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it) 
SilverLife's avatar
SilverLife committed
205
206
207
208
209
            {
                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
210
                ss << std::endl;
SilverLife's avatar
SilverLife committed
211
212
213
214
215
            }
        }

        if (node)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
216
            ss << "Nodes:" << std::endl;
217
            for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) 
SilverLife's avatar
SilverLife committed
218
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
219
220
                ss << rank << ": " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
221
222
223
224
225
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";

                {
226
                    ss << "  " << m->GetMarker(*it,m->NewMarker());
227
228

                    ss << "(" << 
Kirill Terekhov's avatar
Kirill Terekhov committed
229
230
231
                        std::setw(3) << it->RealArray(m->CoordsTag())[0] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[1] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[2] << ")";
232

SilverLife's avatar
SilverLife committed
233
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
234
                ss << std::endl;
SilverLife's avatar
SilverLife committed
235
236
237
            }
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
238
239
        ss << "=========================================" << std::endl;
        os << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
240
    }
241
	*/
SilverLife's avatar
SilverLife committed
242

243
	/*
SilverLife's avatar
SilverLife committed
244
245
246
247
248
    void AdaptiveMesh::UpdateStatus()
    {

        for(ElementType mask = CELL; mask >= NODE; mask = PrevElementType(mask))
        {
249
            for(Mesh::iteratorElement it = m->BeginElement(mask); it != m->EndElement(); it++)
SilverLife's avatar
SilverLife committed
250
251
252
253
254
255
256
257
258
            {
                int stat = 0;
                if (it->GetStatus() == Element::Shared) stat = 1;
                else if (it->GetStatus() == Element::Ghost)  stat = 2;

                tag_status[it->self()] = stat;
            }
        }
    }
259
	*/
SilverLife's avatar
SilverLife committed
260
    
261
262
	
	void AdaptiveMesh::PrintSet(std::ostream & fout, ElementSet set)
SilverLife's avatar
SilverLife committed
263
    {
264
265
266
		fout << "set " << set.GetName();
		fout << " size " << set.Size();
		if(set.HaveParent()) fout << " parent " << set.GetParent().GetName();
267
268
269
		fout << " elements ";
		for(ElementSet::iterator it = set.Begin(); it != set.End(); ++it)
			fout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << ":" << it->GlobalID() << " ";
270
		fout << " children ";
SilverLife's avatar
SilverLife committed
271
		for(ElementSet child = set.GetChild(); child.isValid(); child = child.GetSibling())
272
273
            fout << child.GetName() << " ";
		fout << std::endl;
SilverLife's avatar
SilverLife committed
274
    }
275
	
SilverLife's avatar
SilverLife committed
276

277
	/*
278
279
280
    void AdaptiveMesh::SynchronizeSet(ElementSet set)
    {
    #ifdef USE_MPI
281
282
        int size = m->GetProcessorsNumber();
        int rank = m->GetProcessorRank();
283
284
        for (int i = 0; i < size; i++)
        {
285
286
            set.IntegerArray(m->SendtoTag()).push_back(i);
            m->ExchangeMarked();
287
288
289
        }
    #endif
    }
290
	*/
291

292
	/*
SilverLife's avatar
SilverLife committed
293
294
    void AdaptiveMesh::Test()
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
295
        std::cout << rank << ": ================" << std::endl;
SilverLife's avatar
SilverLife committed
296
297
        PrintSet(root,"");
    }
298
	*/
299
300
301
	
	void AdaptiveMesh::ClearData()
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
302
		set_id        = m->DeleteTag(set_id);
Kirill Terekhov's avatar
Kirill Terekhov committed
303
304
305
		level         = m->DeleteTag(level);
		hanging_nodes = m->DeleteTag(hanging_nodes);
		parent_set    = m->DeleteTag(parent_set);
306
307
308
309
310
311
312
313
		root.DeleteSetTree();
	}
	
	void AdaptiveMesh::PrepareSet()
	{
		//retrive set for coarsening, initialize set if is not present
		if( !root.isValid() )
		{
314
			root = m->GetSet("AM_ROOT_SET");
315
316
			if( root == InvalidElement() )
			{
317
318
				root = m->CreateSetUnique("AM_ROOT_SET").first;
				//root.SetExchange(ElementSet::SYNC_ELEMENTS_SHARED);
319
				level[root] = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
320
				for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
321
322
323
324
				{
					root.PutElement(it->self());
					parent_set[it->self()] = root.GetHandle();
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
325
				m->Enumerate(CELL,set_id);
326
			}
SilverLife's avatar
SilverLife committed
327
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
328
		if( !m->HaveGlobalID(CELL) ) m->AssignGlobalID(CELL); //for unique set names
Kirill Terekhov's avatar
Kirill Terekhov committed
329
		m->ResolveSets();
330
331
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
332
	AdaptiveMesh::AdaptiveMesh(Mesh & _m) : m(&_m)
333
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
334
		model = NULL;
335
		//create a tag that stores maximal refinement level of each element
Kirill Terekhov's avatar
Kirill Terekhov committed
336
		level = m->CreateTag("REFINEMENT_LEVEL",DATA_INTEGER,CELL|FACE|EDGE|NODE|ESET,NONE,1);
337
		//tag_status = m->CreateTag("TAG_STATUS",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
338
		set_id = m->CreateTag("SET_ID",DATA_INTEGER,CELL|ESET,ESET,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
339
340
		//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);
341
		//create a tag that stores links to all the hanging nodes of the cell
Kirill Terekhov's avatar
Kirill Terekhov committed
342
		hanging_nodes = m->CreateTag("HANGING_NODES",DATA_REFERENCE,CELL|FACE,NONE);
343
		//create a tag that stores links to sets
Kirill Terekhov's avatar
Kirill Terekhov committed
344
		parent_set = m->CreateTag("PARENT_SET",DATA_REFERENCE,CELL,NONE,1);
345
346
	    size = m->GetProcessorsNumber();
    	rank = m->GetProcessorRank();
347
	}
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369

	void ReportSub(ElementSet root, int tab, std::fstream & fout)
	{
		if( root.HaveChild() )
		{
			for(ElementSet jt = root.GetChild(); jt.isValid(); jt = jt.GetSibling())
				ReportSub(jt,tab+1,fout);
		}
		Storage::integer_array parr = root.IntegerArray(root.GetMeshLink()->ProcessorsTag());
		for(int q = 0; q < tab; ++q) fout << "  ";
		fout << "set " << root.GetName() << " procs ";
		for(unsigned q = 0; q < parr.size(); ++q) fout << parr[q] << " ";
		fout << " owner " << root.Integer(root.GetMeshLink()->OwnerTag());
		fout << " kids " << root.Size();
		fout << " status " << Element::StatusName(root.GetStatus()) << std::endl;
	}

	void AdaptiveMesh::ReportSets(std::fstream & fout)
	{
		for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it)
			if( !it->HaveParent() ) ReportSub(it->self(),0,fout);
	}
370
371
372
373
374
375
376
	
	AdaptiveMesh::~AdaptiveMesh()
	{
		//do not delete tags, user may want to repetitively use this class
		//as extension of class mesh in limited code span
	}
	
377
	void AdaptiveMesh::CheckParentSet(std::string file, int line)//, TagInteger indicator)
Kirill Terekhov's avatar
Kirill Terekhov committed
378
	{
379
		ENTER_FUNC();
380
#if !defined(NDEBUG)
Kirill Terekhov's avatar
Kirill Terekhov committed
381
		int err = 0;
382
		Storage::integer_array procs;
Kirill Terekhov's avatar
Kirill Terekhov committed
383
384
385
386
		for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
		{
			if( parent_set[*it] == InvalidHandle() )
			{
387
				REPORT_STR(m->GetProcessorRank() << " parent set not valid on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << " " << level[*it]);
388
				std::cout << m->GetProcessorRank() << " parent set not valid on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << " " << level[*it] << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
389
390
391
392
				err++;
			}
			else if( GetHandleElementType(parent_set[*it]) != ESET )
			{
393
				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]);
394
				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;
Kirill Terekhov's avatar
Kirill Terekhov committed
395
396
				err++;
			}
397
			else if( parent_set[*it] != root.GetHandle() )//&& (!indicator.isValid() || indicator[*it]) )
398
399
400
401
402
403
			{
				ElementSet set(m,parent_set[*it]);
				if( !set.HaveParent() )
				{
					REPORT_STR(m->GetProcessorRank() << 
					" parent set " << ElementTypeName(GetHandleElementType(parent_set[*it])) << ":" << GetHandleID(parent_set[*it]) << 
404
					" name " << set.GetName() << " owner " << set.Integer(m->OwnerTag()) << " status " << Element::StatusName(set.GetStatus()) << " kids " << set.Size() <<
405
					" does not have parent " <<
406
					" on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << " owner " << it->Integer(m->OwnerTag()) << " lvl " << level[*it]);
407
408
					std::cout << m->GetProcessorRank() << 
					" parent set " << ElementTypeName(GetHandleElementType(parent_set[*it])) << ":" << GetHandleID(parent_set[*it]) << 
409
					" name " << set.GetName() << " owner " << set.Integer(m->OwnerTag()) << " status " << Element::StatusName(set.GetStatus()) << " kids " << set.Size() <<
410
					" does not have parent " <<
411
412
413
414
415
					" on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << " owner " << it->Integer(m->OwnerTag()) << " lvl " << level[*it] << std::endl;
					procs = it->IntegerArray(m->ProcessorsTag());
					std::cout << "cell procs:"; for(unsigned q = 0; q < procs.size(); ++q) std::cout << " " << procs[q]; std::cout << std::endl;
					procs = set.IntegerArray(m->ProcessorsTag());
					std::cout << "eset procs:"; for(unsigned q = 0; q < procs.size(); ++q) std::cout << " " << procs[q]; std::cout << std::endl;
416
417
418
419
420
421
422
423
424
425
426
427
					err++;
				}
				else
				{
					HandleType parent = set.GetParent().GetHandle();
					if( GetHandleElementType(parent) != ESET )
					{
						REPORT_STR(m->GetProcessorRank() << 
						" parent set " << ElementTypeName(GetHandleElementType(parent_set[*it])) << ":" << GetHandleID(parent_set[*it]) << 
						" name " << set.GetName() << " owner " << set.Integer(m->OwnerTag()) << " status " << Element::StatusName(set.GetStatus()) << 
						" has parent " << ElementTypeName(GetHandleElementType(parent)) << ":" << GetHandleID(parent) <<
						" on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << " " << level[*it]);
428
429
430
431
432
						std::cout << m->GetProcessorRank() << 
						" parent set " << ElementTypeName(GetHandleElementType(parent_set[*it])) << ":" << GetHandleID(parent_set[*it]) << 
						" name " << set.GetName() << " owner " << set.Integer(m->OwnerTag()) << " status " << Element::StatusName(set.GetStatus()) << 
						" has parent " << ElementTypeName(GetHandleElementType(parent)) << ":" << GetHandleID(parent) <<
						" on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << " " << level[*it] << std::endl;
433
434
435
436
						err++;
					}
				}
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
437
438
		}
		err = m->Integrate(err);
439
		EXIT_FUNC();
Kirill Terekhov's avatar
Kirill Terekhov committed
440
441
		if( err ) 
		{
442
443
			REPORT_STR(rank << " error in " << __FUNCTION__ << " " << file << ":" << line);
			std::cout << rank << " error in " << __FUNCTION__ << " " << file << ":" << line << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
444
445
			exit(-1);
		}
446
#endif //NDEBUG
Kirill Terekhov's avatar
Kirill Terekhov committed
447
448
	}
	
449
450
	bool AdaptiveMesh::Refine(TagInteger & indicator)
	{
451
452
		static int fi = 0;
        ENTER_FUNC();
453
454
455
		static int call_counter = 0;
		int ret = 0; //return number of refined cells
		//initialize tree structure
456
457
		//m->CheckCentroids(__FILE__,__LINE__);
		ENTER_BLOCK();
458
		PrepareSet();
459
460
		EXIT_BLOCK();

461
462
		//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
463
		
464
465
466
467
		//ENTER_BLOCK();
		//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
		//	if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
		//EXIT_BLOCK();
468

469
		m->CheckSetLinks(__FILE__,__LINE__);
470
		CheckParentSet(__FILE__,__LINE__);
471
		/*
472
		ENTER_BLOCK();
473
		m->ResolveSets();
474
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
475
		m->ExchangeData(parent_set,CELL,0);
476
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
477
		m->ResolveSets();
478
		//m->CheckCentroids(__FILE__,__LINE__);
479
		//m->ExchangeData(parent_set,CELL,0);
480
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
481
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
482
		//m->CheckCentroids(__FILE__,__LINE__);
483
		//CheckParentSet(Tag());
484
		EXIT_BLOCK();
485
		*/
486
487
488
489
		//ENTER_BLOCK();
		//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
		//	if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
		//EXIT_BLOCK();
490

491
492
		//m->Save("before_refine_parent"+std::to_string(fi)+".pvtk");
		//std::cout << "Save before_refine_parent"+std::to_string(fi)+".pvtk" << std::endl;
493
494
495
		

		
496
497
		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
498
		ENTER_BLOCK();
499
		//0. Extend indicator for edges and faces
Kirill Terekhov's avatar
Kirill Terekhov committed
500
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
501
502
		//1.Communicate indicator - it may be not synced
		m->ExchangeData(indicator,CELL,0);
503
504
		while(scheduled)
		{
505
			REPORT_VAL("scheduled",scheduled);
506
507
			//2.Propogate indicator down to the faces,edges
			//  select schedule for them
508
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
509
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
510
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
511
				Cell c = m->CellByLocalID(it);
512
513
514
515
516
517
518
519
520
521
				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
					}
				}
			}
522
			EXIT_BLOCK();
523
			//3.Communicate indicator on faces and edges
524
			m->ReduceData(indicator,FACE|EDGE,0,ReduceMax);
Kirill Terekhov's avatar
Kirill Terekhov committed
525
			m->ExchangeData(indicator,FACE|EDGE,0);
526
527
528
529
			//4.Check for each cell if there is
			//  any hanging node with adjacent in a need to refine,
			//  schedule for refinement earlier.
			scheduled = 0;
530
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
531
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
532
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
533
				Cell c = m->CellByLocalID(it);
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
				//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;
						}
					}
				}
			}
555
			EXIT_BLOCK();
556
557
558
559
			//5.Exchange indicator on cells
			m->ReduceData(indicator,CELL,0,ReduceMax);
			m->ExchangeData(indicator,CELL,0);
			//6.Go back to 1 until no new elements scheduled
Kirill Terekhov's avatar
Kirill Terekhov committed
560
			scheduled = m->Integrate(scheduled);
561
562
			if( scheduled ) schedule_counter++;
		}
563
		//m->ExchangeData(indicator,CELL | FACE | EDGE,0);
564
		EXIT_BLOCK();
565
566
		
		//m->ExchangeData(hanging_nodes,CELL|FACE,NONE);
Kirill Terekhov's avatar
Kirill Terekhov committed
567
		//m->Save("indicator.pmf");
568
569
570
571
572
573
574
575
		
		//ENTER_BLOCK();
		//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
		//	if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
		//EXIT_BLOCK();
		
		//if( !Element::CheckConnectivity(m) ) std::cout << __FILE__ << ":" << __LINE__ << " broken connectivity" << std::endl;
		//m->CheckCentroids(__FILE__,__LINE__);
576
		//6.Refine
577
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
578
		m->BeginModification();
579
580
581
582
		while(schedule_counter)
		{
			Storage::real xyz[3] = {0,0,0};
			//7.split all edges of the current schedule
583
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
584
			for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
585
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
586
				Edge e = m->EdgeByLocalID(it);
587
588
				if( !e.Hidden() && indicator[e] == schedule_counter )
				{
589
590
591
592
593
594
595
					//ENTER_BLOCK();
					//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
					//	if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
					//EXIT_BLOCK();
					//REPORT_STR("split edge " << e.LocalID() << " nodes " << e.getBeg().LocalID() << "," << e.getEnd().LocalID() << " level " << level[e] << " lc size " << m->LowConn(e.GetHandle()).size() );
					//ElementArray<Node> nodes = e.getNodes();
					//for(int q = 0; q < nodes.size(); ++q) REPORT_STR("node " << nodes[q].GetHandle() << " " << nodes[q].LocalID() << (nodes[q].Hidden()?" hidden " : " good ") );
596
597
598
					//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
599
					for(Storage::integer d = 0; d < m->GetDimensions(); ++d)
600
601
602
						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
603
					Node n = m->CreateNode(xyz);
604
					//set increased level for new node
605
					level[n] = level[e.getBeg()] = level[e.getEnd()] = level[e]+1;
606
607
608
609
					//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
610
					ElementArray<Edge> new_edges = Edge::SplitEdge(e,ElementArray<Node>(m,1,n.GetHandle()),0);
611
612
					//set increased level for new edges
					level[new_edges[0]] = level[new_edges[1]] = level[e]+1;
613
614
615
616
617
618
619
					
					//for(int q = 0; q < 2; ++q)
					//{
					//	REPORT_STR("new edges["<<q<<"]" << new_edges[q].LocalID() << " nodes " << new_edges[q].getBeg().LocalID() << "," << new_edges[q].getEnd().LocalID() << " level " << level[new_edges[q]]);
					//}
					
					//if( !Element::CheckConnectivity(m) ) std::cout << __FILE__ << ":" << __LINE__ << " broken connectivity" << std::endl;
620
621
				}
			}
622
			EXIT_BLOCK();
623
			//8.split all faces of the current schedule, using hanging nodes on edges
624
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
625
			for(Storage::integer it = 0; it < m->FaceLastLocalID(); ++it) if( m->isValidFace(it) )
626
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
627
				Face f = m->FaceByLocalID(it);
628
629
				if( !f.Hidden() && indicator[f] == schedule_counter )
				{
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
#if !defined(NDEBUG)
					ElementArray<Edge> face_edges = f.getEdges();
					for(ElementArray<Edge>::iterator jt = face_edges.begin(); jt != face_edges.end(); ++jt)
					{
						if( level[*jt] != level[f]+1 )
						{
							std::cout << m->GetProcessorRank() << " face " << f.LocalID();
							std::cout << " " << Element::StatusName(f.GetStatus()) << " owner " << f.Integer(m->OwnerTag());
							std::cout << " lvl " << level[f] << " ind " << indicator[f];
							std::cout << " edge " << jt->LocalID();
							std::cout << " " << Element::StatusName(jt->GetStatus()) << " owner " << jt->Integer(m->OwnerTag());
							std::cout << " lvl " << level[*jt] << " ind " << indicator[*jt];
							std::cout << std::endl;
							ElementArray<Cell> adj_cells;
							adj_cells = f.getCells();
							std::cout << "face cells ";
							for(ElementArray<Cell>::iterator kt = adj_cells.begin(); kt != adj_cells.end(); ++kt)
							{
								std::cout << m->GetProcessorRank() << " face " << kt->LocalID();
								std::cout << " " << Element::StatusName(kt->GetStatus()) << " owner " << kt->Integer(m->OwnerTag());
								std::cout << " lvl " << level[*kt] << " ind " << indicator[*kt];
								std::cout << std::endl;
							}
							std::cout << "edge cells ";
							adj_cells = jt->getCells();
							for(ElementArray<Cell>::iterator kt = adj_cells.begin(); kt != adj_cells.end(); ++kt)
							{
								std::cout << m->GetProcessorRank() << " face " << kt->LocalID();
								std::cout << " " << Element::StatusName(kt->GetStatus()) << " owner " << kt->Integer(m->OwnerTag());
								std::cout << " lvl " << level[*kt] << " ind " << indicator[*kt];
								std::cout << std::endl;
							}
						}
						assert(level[*jt] == level[f]+1);
					}
#endif //NDEBUG
666
667
668
669
670
671
672
673
674
675
676
677
678
					//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
679
					Node n = m->CreateNode(xyz);
680
					//set increased level for the new node
681
					level[n] = level[f]+1;
682
683
684
					//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
685
686
					ElementArray<Node> edge_nodes(m,2); //to create new edges
					ElementArray<Edge> hanging_edges(m,face_hanging_nodes.size());
687
688
689
690
					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
691
						hanging_edges[kt] = m->CreateEdge(edge_nodes).first;
692
693
694
695
696
697
698
699
700
701
						//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;
				}
			}
702
			EXIT_BLOCK();
703
			//this tag helps recreate internal face
Kirill Terekhov's avatar
Kirill Terekhov committed
704
			TagReferenceArray internal_face_edges = m->CreateTag("INTERNAL_FACE_EDGES",DATA_REFERENCE,NODE,NODE,4);
705
			//this marker helps detect edges of current cell only
Kirill Terekhov's avatar
Kirill Terekhov committed
706
			MarkerType mark_cell_edges = m->CreateMarker();
707
			//this marker helps detect nodes hanging on edges of unrefined cell
Kirill Terekhov's avatar
Kirill Terekhov committed
708
			MarkerType mark_hanging_nodes = m->CreateMarker();
709
			//9.split all cells of the current schedule
710
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
711
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
712
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
713
				Cell c = m->CellByLocalID(it);
714
715
716
717
718
719
720
721
722
723
724
				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
725
					Node n = m->CreateNode(xyz);
726
					//set increased level for the new node
727
					level[n] = level[c]+1;
728
729
					//retrive all edges of current face to mark them
					ElementArray<Edge> cell_edges = c.getEdges();
730
731
732
733
734
#if !defined(NDEBUG)
					for(ElementArray<Edge>::iterator jt = cell_edges.begin(); jt != cell_edges.end(); ++jt) assert(level[*jt] == level[c]+1);
					ElementArray<Face> cell_faces = c.getFaces();
					for(ElementArray<Face>::iterator jt = cell_faces.begin(); jt != cell_faces.end(); ++jt) assert(level[*jt] == level[c]+1);
#endif //NDEBUG
735
736
737
					//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
738
739
					ElementArray<Node> edge_nodes(m,2);
					ElementArray<Edge> edges_to_faces(m,cell_hanging_nodes.size());
740
741
742
743
744
745
					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
746
						edges_to_faces[kt] = m->CreateEdge(edge_nodes).first;
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
						//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
782
					ElementArray<Face> internal_faces(m,edge_hanging_nodes.size());
783
784
785
786
787
788
789
790
791
792
					//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
793
						internal_faces[kt] = m->CreateFace(ElementArray<Edge>(m,face_edges.begin(),face_edges.end())).first;
794
795
796
797
798
						//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
799
800
801
802
803
804
805
					//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;
					//	
					//}
806
807
					//split the cell
					//retrive parent set
Kirill Terekhov's avatar
Kirill Terekhov committed
808
					ElementSet parent(m,parent_set[c]);
809
					//create set corresponding to old coarse cell
810
811
					//Storage::real cnt[3];
					//c.Centroid(cnt);
812
					std::stringstream set_name;
Kirill Terekhov's avatar
Kirill Terekhov committed
813
					//set_name << parent.GetName() << "_C" << c.GlobalID(); //rand may be unsafe
Kirill Terekhov's avatar
Kirill Terekhov committed
814
					if( parent == root )
815
						set_name << "AM_R" << set_id[c];
Kirill Terekhov's avatar
Kirill Terekhov committed
816
					else
817
						set_name << parent.GetName() << "C" << set_id[c];
818
					//set_name << base64_encode_((unsigned char *)cnt,3*sizeof(double)/sizeof(unsigned char));
Kirill Terekhov's avatar
Kirill Terekhov committed
819
					
Kirill Terekhov's avatar
Kirill Terekhov committed
820
821
822
823
824
825
					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;
826
						std::cout << rank << " Elements of " << check_set.GetName() << ": ";
Kirill Terekhov's avatar
Kirill Terekhov committed
827
						for(ElementSet::iterator it = check_set.Begin(); it != check_set.End(); ++it)
828
							std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << "," << indicator[*it] << " ";
Kirill Terekhov's avatar
Kirill Terekhov committed
829
						std::cout << std::endl;
830
831
832
833
834
835
836
837
838
839
840
						std::cout << rank << " Elements of " << parent.GetName() << ": ";
						for(ElementSet::iterator it = parent.Begin(); it != parent.End(); ++it)
							std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << "," << indicator[*it] << " ";
						std::cout << std::endl;
						if( parent.HaveChild() )
						{
							std::cout << rank << " Children of " << parent.GetName() << ": ";
							for(ElementSet jt = parent.GetChild(); jt.isValid(); jt = jt.GetSibling() )
								std::cout << jt.GetName() << " size " << jt.Size() << " ";
							std::cout << std::endl;
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
841
						exit(-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
842
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
843
					
Kirill Terekhov's avatar
Kirill Terekhov committed
844
					ElementSet cell_set = m->CreateSetUnique(set_name.str()).first;
845
					//cell_set->SetExchange(ElementSet::SYNC_ELEMENTS_ALL);
846
					level[cell_set] = level[c]+1;
847
848
849
850
851
					set_id[cell_set] = set_id[c];

					ElementArray<Cell> new_cells = Cell::SplitCell(c,internal_faces,0);
					std::sort(new_cells.begin(),new_cells.end(),Mesh::CentroidComparator(m));
					
852
853
854
					//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
855
						set_id[new_cells[kt]] = kt;
856
857
858
859
						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
860
861
862
863
864
865
866
867
868
869
					/*
					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() )
870
					parent.AddChild(cell_set);
Kirill Terekhov's avatar
Kirill Terekhov committed
871
					//else assert(cell_set->GetParent() == parent);
872
873
874
875
					//increment number of refined cells
					ret++;
				}
			}
876
			EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
877
878
879
			m->ReleaseMarker(mark_hanging_nodes);
			m->ReleaseMarker(mark_cell_edges);
			m->DeleteTag(internal_face_edges);
880
881
882
			//10.jump to later schedule, and go to 7.
			schedule_counter--;
		}
883
		m->CheckSetLinks(__FILE__,__LINE__);
884
		//free created tag
Kirill Terekhov's avatar
Kirill Terekhov committed
885
		m->DeleteTag(indicator,FACE|EDGE);
SilverLife's avatar
SilverLife committed
886
887


888
		m->CheckSetLinks(__FILE__,__LINE__);
889
		//11. Restore parallel connectivity, global ids
Kirill Terekhov's avatar
Kirill Terekhov committed
890
		m->ResolveModification();
Kirill Terekhov's avatar
Kirill Terekhov committed
891
		//m->SynchronizeMarker(m->NewMarker(),CELL|FACE|EDGE|NODE,SYNC_BIT_OR);
892
        //ExchangeGhost(3,NODE); // Construct Ghost cells in 2 layers connected via nodes
893
		//12. Let the user update their data
894
		//todo: call back function for New( cells
Kirill Terekhov's avatar
Kirill Terekhov committed
895
		if( model ) model->Adaptation(*m);
896
		m->CheckSetLinks(__FILE__,__LINE__);
897
		//13. Delete old elements of the mesh
Kirill Terekhov's avatar
Kirill Terekhov committed
898
		m->ApplyModification();
Kirill Terekhov's avatar
Kirill Terekhov committed
899
		
Kirill Terekhov's avatar
Kirill Terekhov committed
900
		//m->ExchangeGhost(1,NODE,m->NewMarker());
901
		//14. Done
Kirill Terekhov's avatar
Kirill Terekhov committed
902
        //cout << rank << ": Before end " << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
903
		m->EndModification();
904
		EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
905
		
906
		m->ExchangeData(parent_set,CELL,0);
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931

		/*
		ENTER_BLOCK();
		for(Storage::integer it = 0; it < m->EsetLastLocalID(); ++it) if( m->isValidElementSet(it) )
		{
			ElementSet set = m->EsetByLocalID(it);
			if( set.GetName().substr(0,3) == "AM_" )
			{
				//int imax = -1, imin = INT_MAX;
				//for(ElementSet::iterator jt = set.Begin(); jt != set.End(); ++jt)
				//{
				//	imax = std::max(imax,indicator[*jt]);
				//	imin = std::min(imin,indicator[*jt]);
				//}
				//std::cout << "on proc " << m->GetProcessorRank() << " set " << set.GetName() << " size " << set.Size() << " set indicator " << indicator[set] << " elements indicator " << imin << ":" << imax;
				//if( set.HaveParent() ) std::cout << " parent " << set.GetParent().GetName();
				//std::cout << std::endl;
				//if( !set.HaveChild() )
					set.SynchronizeSetParents();
			}
		}
		m->ExchangeMarked();
		EXIT_BLOCK();
		*/
		CheckParentSet(__FILE__,__LINE__);
932
		
933
934
		//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
935
		fi++;
Kirill Terekhov's avatar
Kirill Terekhov committed
936
		
937
		//ExchangeData(hanging_nodes,CELL | FACE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
938
        //m->ResolveSets();
939

Kirill Terekhov's avatar
Kirill Terekhov committed
940
941
942
943
944
945
        //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
946
    	//m->ResolveSets();
Kirill Terekhov's avatar
Kirill Terekhov committed
947
		
Kirill Terekhov's avatar
Kirill Terekhov committed
948
		//restore face orientation
Kirill Terekhov's avatar
Kirill Terekhov committed
949
		//BUG: bad orientation not fixed automatically
950
951
		
		/*
952
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
953
954
955
956
957
958
959
960
		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;
961
962
		if( nfixed ) REPORT_STR(rank << " fixed " << nfixed << " faces");
		EXIT_BLOCK();
963
964
		 */
		
Kirill Terekhov's avatar
Kirill Terekhov committed
965
966
967
968
969
970
971
972
973
		/*
		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
974
975
		//std::cout << rank << " check parent_set at end" << std::endl;
		//CheckParentSet();
976
977

        //ExchangeData(hanging_nodes,CELL | FACE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
978
        //cout << rank << ": After end " << std::endl;
979
		//reorder element's data to free up space