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

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

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

Kirill Terekhov's avatar
Kirill Terekhov committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
	std::string base64_encode_(unsigned char const* buf, unsigned int bufLen)
	{
		static const std::string base64_chars =
		"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
		"abcdefghijklmnopqrstuvwxyz"
		"0123456789+/";
		std::string ret;
		int i = 0;
		int j = 0;
		unsigned char char_array_3[3];
		unsigned char char_array_4[4];
		
		while (bufLen--)
		{
			char_array_3[i++] = *(buf++);
			if (i == 3)
			{
				char_array_4[0] = (char_array_3[0] & 0xfc) >> 2;
				char_array_4[1] = ((char_array_3[0] & 0x03) << 4) + ((char_array_3[1] & 0xf0) >> 4);
				char_array_4[2] = ((char_array_3[1] & 0x0f) << 2) + ((char_array_3[2] & 0xc0) >> 6);
				char_array_4[3] = char_array_3[2] & 0x3f;
				for(i = 0; (i <4) ; i++)
					ret += base64_chars[char_array_4[i]];
				i = 0;
			}
		}
		
		if (i)
		{
			for(j = i; j < 3; j++)
				char_array_3[j] = '\0';
	  
			char_array_4[0] = (char_array_3[0] & 0xfc) >> 2;
			char_array_4[1] = ((char_array_3[0] & 0x03) << 4) + ((char_array_3[1] & 0xf0) >> 4);
			char_array_4[2] = ((char_array_3[1] & 0x0f) << 2) + ((char_array_3[2] & 0xc0) >> 6);
			char_array_4[3] = char_array_3[2] & 0x3f;
	  
			for (j = 0; (j < i + 1); j++)
				ret += base64_chars[char_array_4[j]];
	  
			while((i++ < 3))
				ret += '=';
		}
		
		return ret;
	}
87

88
89
90
/// todo:
/// 1. coarsment
/// 2. strategy for faces/edges with faults
Kirill Terekhov's avatar
Kirill Terekhov committed
91
/// 3. geom model supportbn
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
/// 4. make class abstract virtual for user implementation of refinement and coarsment indicators
/// see in code todo:
namespace INMOST
{
	void CleanupSets(ElementSet set)
	{
		ElementSet::iterator it = set.Begin();
		while(it != set.End())
		{
			if( it->isValid() ) ++it;
			else it = set.Erase(it);
		}
		for(ElementSet child = set.GetChild(); child.isValid(); child = child.GetSibling())
			CleanupSets(child);
	}
	
	void ReduceMax(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void) size;
		element->Integer(tag) = std::max(element->Integer(tag),*((const INMOST_DATA_INTEGER_TYPE *)data));
	}
	
	void ReduceMin(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void) size;
		element->Integer(tag) = std::min(element->Integer(tag),*((const INMOST_DATA_INTEGER_TYPE *)data));
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
119
120
121
122
123
124
125
126
127
128
	
	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
129

130
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
131
    void AdaptiveMesh::PrintSetLocal(std::string offset, ElementSet it, std::stringstream& ss)
SilverLife's avatar
SilverLife committed
132
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
133
134
        std::stringstream ss1;
        ss1 << offset << rank << ": Set : " << std::setw(5) << it.GetName() << " ";
SilverLife's avatar
SilverLife committed
135
136
        for (int i = ss1.str().length(); i < 23; i++) ss1 << " ";
        ss << ss1.str();
Kirill Terekhov's avatar
Kirill Terekhov committed
137
        ss << std::setw(6);
SilverLife's avatar
SilverLife committed
138
139
140
141
142
        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";

143
        ss << " tag_owner (" << it.IntegerDF(m->OwnerTag()) << ")";
SilverLife's avatar
SilverLife committed
144
145
146

        //ss << "   level (" << level[it.self()] << ")  ";
        ss << " tag_processors (";
Kirill Terekhov's avatar
Kirill Terekhov committed
147
        std::stringstream ss2;
148
        Storage::integer_array arr = it.IntegerArrayDV(m->ProcessorsTag());
SilverLife's avatar
SilverLife committed
149
150
        for (int i = 0; i < arr.size(); i++)
            ss2 << arr[i] << " ";
Kirill Terekhov's avatar
Kirill Terekhov committed
151
        ss << std::setw(5) << ss2.str() <<")";
SilverLife's avatar
SilverLife committed
152
153
154
155
156
157
158
159
        

        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
160
            std::string type = "unknw";
SilverLife's avatar
SilverLife committed
161
162
163
164
            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
165
            ss << type << "-" << std::setw(2) << p->GlobalID() << " ";
SilverLife's avatar
SilverLife committed
166
167
            p++;
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
168
        ss << std::endl;
SilverLife's avatar
SilverLife committed
169
170
171
172
173
174

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

177
	/*
SilverLife's avatar
SilverLife committed
178
179
    void AdaptiveMesh::PrintSet()
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
180
        std::stringstream ss;
181
        for(Mesh::iteratorSet it = m->BeginSet(); it != m->EndSet(); ++it) 
SilverLife's avatar
SilverLife committed
182
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
183
            //if (it->HaveParent()) continue;
184
            PrintSetLocal("",ElementSet(m,*it),ss);
SilverLife's avatar
SilverLife committed
185
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
186
        std::cout << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
187
    }
188
	*/
SilverLife's avatar
SilverLife committed
189

190
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
191
    void PrintRefs(std::ostream& os, Storage::reference_array refs)
192
193
194
    {
        for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
195
            std::string type = "unknw";
196
197
198
199
200
201
202
            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
203
            os << "(" << xyz[0] << "," << xyz[1] << "," << xyz[2] <<")" <<  std::endl;
204
205
        }
    }
206
	*/
207

208
	/*
Kirill Terekhov's avatar
Kirill Terekhov committed
209
    void AdaptiveMesh::PrintMesh(std::ostream& os, int cell, int face, int edge, int node)
SilverLife's avatar
SilverLife committed
210
211
    {
        if (cell + face + edge + node == 0) return;
Kirill Terekhov's avatar
Kirill Terekhov committed
212
213
        std::stringstream ss;
        ss << "================= " << rank << " =====================" << std::endl;
SilverLife's avatar
SilverLife committed
214
215
        if (cell)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
216
            ss << "Cells: " << m->NumberOfCells() <<  std::endl;
217
            for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) 
SilverLife's avatar
SilverLife committed
218
            {
219
                ss << rank << ": " << it->GlobalID() << " - " << it->LocalID() << " - ";            
SilverLife's avatar
SilverLife committed
220
221
222
                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
223
                ss << std::endl;
SilverLife's avatar
SilverLife committed
224
225
226
227
228
            }
        }

        if (face)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
229
            ss << "Faces: " << m->NumberOfFaces() <<  std::endl;
230
            for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) 
SilverLife's avatar
SilverLife committed
231
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
232
233
                ss << rank << ": " << std::setw(2) << it->LocalID() << " " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
234
235
236
237
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";

238
239
240

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

243
                ss << "  " << m->GetMarker(*it,m->NewMarker());
244
245
246
247

                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
248
                    ss << std::setw(2) << node->GlobalID() << " ";
249
                ss << ")";
Kirill Terekhov's avatar
Kirill Terekhov committed
250
                ss << std::endl;
SilverLife's avatar
SilverLife committed
251
252
253
254
255
            }
        }

        if (edge)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
256
            ss << "Edges: " << m->NumberOfEdges() <<  std::endl;
257
            for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it) 
SilverLife's avatar
SilverLife committed
258
259
260
261
262
            {
                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
263
                ss << std::endl;
SilverLife's avatar
SilverLife committed
264
265
266
267
268
            }
        }

        if (node)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
269
            ss << "Nodes:" << std::endl;
270
            for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) 
SilverLife's avatar
SilverLife committed
271
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
272
273
                ss << rank << ": " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
274
275
276
277
278
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";

                {
279
                    ss << "  " << m->GetMarker(*it,m->NewMarker());
280
281

                    ss << "(" << 
Kirill Terekhov's avatar
Kirill Terekhov committed
282
283
284
                        std::setw(3) << it->RealArray(m->CoordsTag())[0] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[1] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[2] << ")";
285

SilverLife's avatar
SilverLife committed
286
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
287
                ss << std::endl;
SilverLife's avatar
SilverLife committed
288
289
290
            }
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
291
292
        ss << "=========================================" << std::endl;
        os << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
293
    }
294
	*/
SilverLife's avatar
SilverLife committed
295

296
	/*
SilverLife's avatar
SilverLife committed
297
298
299
300
301
    void AdaptiveMesh::UpdateStatus()
    {

        for(ElementType mask = CELL; mask >= NODE; mask = PrevElementType(mask))
        {
302
            for(Mesh::iteratorElement it = m->BeginElement(mask); it != m->EndElement(); it++)
SilverLife's avatar
SilverLife committed
303
304
305
306
307
308
309
310
311
            {
                int stat = 0;
                if (it->GetStatus() == Element::Shared) stat = 1;
                else if (it->GetStatus() == Element::Ghost)  stat = 2;

                tag_status[it->self()] = stat;
            }
        }
    }
312
	*/
SilverLife's avatar
SilverLife committed
313
    
314
315
	
	void AdaptiveMesh::PrintSet(std::ostream & fout, ElementSet set)
SilverLife's avatar
SilverLife committed
316
    {
317
318
319
		fout << "set " << set.GetName();
		fout << " size " << set.Size();
		if(set.HaveParent()) fout << " parent " << set.GetParent().GetName();
320
321
322
		fout << " elements ";
		for(ElementSet::iterator it = set.Begin(); it != set.End(); ++it)
			fout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << ":" << it->GlobalID() << " ";
323
		fout << " children ";
SilverLife's avatar
SilverLife committed
324
		for(ElementSet child = set.GetChild(); child.isValid(); child = child.GetSibling())
325
326
            fout << child.GetName() << " ";
		fout << std::endl;
SilverLife's avatar
SilverLife committed
327
    }
328
	
SilverLife's avatar
SilverLife committed
329

330
	/*
331
332
333
    void AdaptiveMesh::SynchronizeSet(ElementSet set)
    {
    #ifdef USE_MPI
334
335
        int size = m->GetProcessorsNumber();
        int rank = m->GetProcessorRank();
336
337
        for (int i = 0; i < size; i++)
        {
338
339
            set.IntegerArray(m->SendtoTag()).push_back(i);
            m->ExchangeMarked();
340
341
342
        }
    #endif
    }
343
	*/
344

345
	/*
SilverLife's avatar
SilverLife committed
346
347
    void AdaptiveMesh::Test()
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
348
        std::cout << rank << ": ================" << std::endl;
SilverLife's avatar
SilverLife committed
349
350
        PrintSet(root,"");
    }
351
	*/
352
353
354
	
	void AdaptiveMesh::ClearData()
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
355
		set_id        = m->DeleteTag(set_id);
Kirill Terekhov's avatar
Kirill Terekhov committed
356
357
358
		level         = m->DeleteTag(level);
		hanging_nodes = m->DeleteTag(hanging_nodes);
		parent_set    = m->DeleteTag(parent_set);
359
360
361
362
363
364
365
366
		root.DeleteSetTree();
	}
	
	void AdaptiveMesh::PrepareSet()
	{
		//retrive set for coarsening, initialize set if is not present
		if( !root.isValid() )
		{
367
			root = m->GetSet("AM_ROOT_SET");
368
369
			if( root == InvalidElement() )
			{
370
371
				root = m->CreateSetUnique("AM_ROOT_SET").first;
				//root.SetExchange(ElementSet::SYNC_ELEMENTS_SHARED);
372
				level[root] = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
373
				for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
374
375
376
377
				{
					root.PutElement(it->self());
					parent_set[it->self()] = root.GetHandle();
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
378
				m->Enumerate(CELL,set_id);
379
			}
SilverLife's avatar
SilverLife committed
380
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
381
		if( !m->HaveGlobalID(CELL) ) m->AssignGlobalID(CELL); //for unique set names
Kirill Terekhov's avatar
Kirill Terekhov committed
382
		m->ResolveSets();
383
384
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
385
	AdaptiveMesh::AdaptiveMesh(Mesh & _m) : m(&_m)
386
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
387
		model = NULL;
388
		//create a tag that stores maximal refinement level of each element
Kirill Terekhov's avatar
Kirill Terekhov committed
389
		level = m->CreateTag("REFINEMENT_LEVEL",DATA_INTEGER,CELL|FACE|EDGE|NODE|ESET,NONE,1);
390
		//tag_status = m->CreateTag("TAG_STATUS",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
391
		set_id = m->CreateTag("SET_ID",DATA_INTEGER,CELL|ESET,ESET,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
392
393
		//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);
394
		//create a tag that stores links to all the hanging nodes of the cell
Kirill Terekhov's avatar
Kirill Terekhov committed
395
		hanging_nodes = m->CreateTag("HANGING_NODES",DATA_REFERENCE,CELL|FACE,NONE);
396
		//create a tag that stores links to sets
Kirill Terekhov's avatar
Kirill Terekhov committed
397
		parent_set = m->CreateTag("PARENT_SET",DATA_REFERENCE,CELL,NONE,1);
398
399
	    size = m->GetProcessorsNumber();
    	rank = m->GetProcessorRank();
400
	}
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422

	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);
	}
423
424
425
426
427
428
429
	
	AdaptiveMesh::~AdaptiveMesh()
	{
		//do not delete tags, user may want to repetitively use this class
		//as extension of class mesh in limited code span
	}
	
430
	void AdaptiveMesh::CheckParentSet(std::string file, int line)//, TagInteger indicator)
Kirill Terekhov's avatar
Kirill Terekhov committed
431
	{
432
		ENTER_FUNC();
Kirill Terekhov's avatar
Kirill Terekhov committed
433
		int err = 0;
434
		Storage::integer_array procs;
Kirill Terekhov's avatar
Kirill Terekhov committed
435
436
437
438
		for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
		{
			if( parent_set[*it] == InvalidHandle() )
			{
439
				REPORT_STR(m->GetProcessorRank() << " parent set not valid on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << " " << level[*it]);
440
				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
441
442
443
444
				err++;
			}
			else if( GetHandleElementType(parent_set[*it]) != ESET )
			{
445
				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]);
446
				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
447
448
				err++;
			}
449
			else if( parent_set[*it] != root.GetHandle() )//&& (!indicator.isValid() || indicator[*it]) )
450
451
452
453
454
455
			{
				ElementSet set(m,parent_set[*it]);
				if( !set.HaveParent() )
				{
					REPORT_STR(m->GetProcessorRank() << 
					" parent set " << ElementTypeName(GetHandleElementType(parent_set[*it])) << ":" << GetHandleID(parent_set[*it]) << 
456
					" name " << set.GetName() << " owner " << set.Integer(m->OwnerTag()) << " status " << Element::StatusName(set.GetStatus()) << " kids " << set.Size() <<
457
					" does not have parent " <<
458
					" on CELL:" << it->LocalID() << " " << Element::StatusName(it->GetStatus()) << " owner " << it->Integer(m->OwnerTag()) << " lvl " << level[*it]);
459
460
					std::cout << m->GetProcessorRank() << 
					" parent set " << ElementTypeName(GetHandleElementType(parent_set[*it])) << ":" << GetHandleID(parent_set[*it]) << 
461
					" name " << set.GetName() << " owner " << set.Integer(m->OwnerTag()) << " status " << Element::StatusName(set.GetStatus()) << " kids " << set.Size() <<
462
					" does not have parent " <<
463
464
465
466
467
					" 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;
468
469
470
471
472
473
474
475
476
477
478
479
					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]);
480
481
482
483
484
						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;
485
486
487
488
						err++;
					}
				}
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
489
490
		}
		err = m->Integrate(err);
491
		EXIT_FUNC();
Kirill Terekhov's avatar
Kirill Terekhov committed
492
493
		if( err ) 
		{
494
495
			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
496
497
498
499
500
			exit(-1);
		}
		
	}
	
501
502
	bool AdaptiveMesh::Refine(TagInteger & indicator)
	{
503
504
		static int fi = 0;
        ENTER_FUNC();
505
506
507
		static int call_counter = 0;
		int ret = 0; //return number of refined cells
		//initialize tree structure
508
509
		//m->CheckCentroids(__FILE__,__LINE__);
		ENTER_BLOCK();
510
		PrepareSet();
511
512
		EXIT_BLOCK();

513
514
		//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
515
		
516
517
518
519
		//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();
520

521
		m->CheckSetLinks(__FILE__,__LINE__);
522
		CheckParentSet(__FILE__,__LINE__);
523
		/*
524
		ENTER_BLOCK();
525
		m->ResolveSets();
526
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
527
		m->ExchangeData(parent_set,CELL,0);
528
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
529
		m->ResolveSets();
530
		//m->CheckCentroids(__FILE__,__LINE__);
531
		//m->ExchangeData(parent_set,CELL,0);
532
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
533
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
534
		//m->CheckCentroids(__FILE__,__LINE__);
535
		//CheckParentSet(Tag());
536
		EXIT_BLOCK();
537
		*/
538
539
540
541
		//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();
542

543
544
		//m->Save("before_refine_parent"+std::to_string(fi)+".pvtk");
		//std::cout << "Save before_refine_parent"+std::to_string(fi)+".pvtk" << std::endl;
545
546
547
		

		
548
549
		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
550
		ENTER_BLOCK();
551
		//0. Extend indicator for edges and faces
Kirill Terekhov's avatar
Kirill Terekhov committed
552
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
553
554
		while(scheduled)
		{
555
			REPORT_VAL("scheduled",scheduled);
556
			//1.Communicate indicator - it may be not synced
Kirill Terekhov's avatar
Kirill Terekhov committed
557
			m->ExchangeData(indicator,CELL,0);
558
559
			//2.Propogate indicator down to the faces,edges
			//  select schedule for them
560
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
561
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
562
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
563
				Cell c = m->CellByLocalID(it);
564
565
566
567
568
569
570
571
572
573
				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
					}
				}
			}
574
			EXIT_BLOCK();
575
			//3.Communicate indicator on faces and edges
Kirill Terekhov's avatar
Kirill Terekhov committed
576
			m->ExchangeData(indicator,FACE|EDGE,0);
577
578
579
580
			//4.Check for each cell if there is
			//  any hanging node with adjacent in a need to refine,
			//  schedule for refinement earlier.
			scheduled = 0;
581
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
582
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
583
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
584
				Cell c = m->CellByLocalID(it);
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
				//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;
						}
					}
				}
			}
606
			EXIT_BLOCK();
607
			//5.Go back to 1 until no new elements scheduled
Kirill Terekhov's avatar
Kirill Terekhov committed
608
			scheduled = m->Integrate(scheduled);
609
610
			if( scheduled ) schedule_counter++;
		}
611
		m->ExchangeData(indicator,CELL | FACE | EDGE,0);
612
		EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
613
		//m->Save("indicator.pmf");
614
615
616
617
618
619
620
621
		
		//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__);
622
		//6.Refine
623
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
624
		m->BeginModification();
625
626
627
628
		while(schedule_counter)
		{
			Storage::real xyz[3] = {0,0,0};
			//7.split all edges of the current schedule
629
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
630
			for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
631
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
632
				Edge e = m->EdgeByLocalID(it);
633
634
				if( !e.Hidden() && indicator[e] == schedule_counter )
				{
635
636
637
638
639
640
641
					//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 ") );
642
643
644
					//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
645
					for(Storage::integer d = 0; d < m->GetDimensions(); ++d)
646
647
648
						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
649
					Node n = m->CreateNode(xyz);
650
					//set increased level for new node
651
					level[n] = level[e.getBeg()] = level[e.getEnd()] = level[e]+1;
652
653
654
655
					//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
656
					ElementArray<Edge> new_edges = Edge::SplitEdge(e,ElementArray<Node>(m,1,n.GetHandle()),0);
657
658
					//set increased level for new edges
					level[new_edges[0]] = level[new_edges[1]] = level[e]+1;
659
660
661
662
663
664
665
					
					//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;
666
667
				}
			}
668
			EXIT_BLOCK();
669
			//8.split all faces of the current schedule, using hanging nodes on edges
670
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
671
			for(Storage::integer it = 0; it < m->FaceLastLocalID(); ++it) if( m->isValidFace(it) )
672
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
673
				Face f = m->FaceByLocalID(it);
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
				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
689
					Node n = m->CreateNode(xyz);
690
					//set increased level for the new node
691
					level[n] = level[f]+1;
692
693
694
					//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
695
696
					ElementArray<Node> edge_nodes(m,2); //to create new edges
					ElementArray<Edge> hanging_edges(m,face_hanging_nodes.size());
697
698
699
700
					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
701
						hanging_edges[kt] = m->CreateEdge(edge_nodes).first;
702
703
704
705
706
707
708
709
710
711
						//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;
				}
			}
712
			EXIT_BLOCK();
713
			//this tag helps recreate internal face
Kirill Terekhov's avatar
Kirill Terekhov committed
714
			TagReferenceArray internal_face_edges = m->CreateTag("INTERNAL_FACE_EDGES",DATA_REFERENCE,NODE,NODE,4);
715
			//this marker helps detect edges of current cell only
Kirill Terekhov's avatar
Kirill Terekhov committed
716
			MarkerType mark_cell_edges = m->CreateMarker();
717
			//this marker helps detect nodes hanging on edges of unrefined cell
Kirill Terekhov's avatar
Kirill Terekhov committed
718
			MarkerType mark_hanging_nodes = m->CreateMarker();
719
			//9.split all cells of the current schedule
720
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
721
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
722
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
723
				Cell c = m->CellByLocalID(it);
724
725
726
727
728
729
730
731
732
733
734
				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
735
					Node n = m->CreateNode(xyz);
736
					//set increased level for the new node
737
					level[n] = level[c]+1;
738
739
740
741
742
					//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
743
744
					ElementArray<Node> edge_nodes(m,2);
					ElementArray<Edge> edges_to_faces(m,cell_hanging_nodes.size());
745
746
747
748
749
750
					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
751
						edges_to_faces[kt] = m->CreateEdge(edge_nodes).first;
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
782
783
784
785
786
						//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
787
					ElementArray<Face> internal_faces(m,edge_hanging_nodes.size());
788
789
790
791
792
793
794
795
796
797
					//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
798
						internal_faces[kt] = m->CreateFace(ElementArray<Edge>(m,face_edges.begin(),face_edges.end())).first;
799
800
801
802
803
						//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
804
805
806
807
808
809
810
					//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;
					//	
					//}
811
812
					//split the cell
					//retrive parent set
Kirill Terekhov's avatar
Kirill Terekhov committed
813
					ElementSet parent(m,parent_set[c]);
814
					//create set corresponding to old coarse cell
815
816
					//Storage::real cnt[3];
					//c.Centroid(cnt);
817
					std::stringstream set_name;
Kirill Terekhov's avatar
Kirill Terekhov committed
818
					//set_name << parent.GetName() << "_C" << c.GlobalID(); //rand may be unsafe
Kirill Terekhov's avatar
Kirill Terekhov committed
819
					if( parent == root )
820
						set_name << "AM_R" << set_id[c];
Kirill Terekhov's avatar
Kirill Terekhov committed
821
					else
822
						set_name << parent.GetName() << "C" << set_id[c];
823
					//set_name << base64_encode_((unsigned char *)cnt,3*sizeof(double)/sizeof(unsigned char));
Kirill Terekhov's avatar
Kirill Terekhov committed
824
					
Kirill Terekhov's avatar
Kirill Terekhov committed
825
826
827
828
829
830
					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;
831
						std::cout << rank << " Elements of " << check_set.GetName() << ": ";
Kirill Terekhov's avatar
Kirill Terekhov committed
832
						for(ElementSet::iterator it = check_set.Begin(); it != check_set.End(); ++it)
833
							std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << "," << indicator[*it] << " ";
Kirill Terekhov's avatar
Kirill Terekhov committed
834
						std::cout << std::endl;
835
836
837
838
839
840
841
842
843
844
845
						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
846
						exit(-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
847
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
848
					
Kirill Terekhov's avatar
Kirill Terekhov committed
849
					ElementSet cell_set = m->CreateSetUnique(set_name.str()).first;
850
					//cell_set->SetExchange(ElementSet::SYNC_ELEMENTS_ALL);
851
					level[cell_set] = level[c]+1;
852
853
854
855
856
					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));
					
857
858
859
					//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
860
						set_id[new_cells[kt]] = kt;
861
862
863
864
						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
865
866
867
868
869
870
871
872
873
874
					/*
					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() )
875
					parent.AddChild(cell_set);
Kirill Terekhov's avatar
Kirill Terekhov committed
876
					//else assert(cell_set->GetParent() == parent);
877
878
879
880
					//increment number of refined cells
					ret++;
				}
			}
881
			EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
882
883
884
			m->ReleaseMarker(mark_hanging_nodes);
			m->ReleaseMarker(mark_cell_edges);
			m->DeleteTag(internal_face_edges);
885
886
887
			//10.jump to later schedule, and go to 7.
			schedule_counter--;
		}
888
		m->CheckSetLinks(__FILE__,__LINE__);
889
		//free created tag
Kirill Terekhov's avatar
Kirill Terekhov committed
890
		m->DeleteTag(indicator,FACE|EDGE);
SilverLife's avatar
SilverLife committed
891
892


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

		/*
		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__);
937
		
938
939
		//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
940
		fi++;
Kirill Terekhov's avatar
Kirill Terekhov committed
941
		
942
		//ExchangeData(hanging_nodes,CELL | FACE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
943
        //m->ResolveSets();
944

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

        //ExchangeData(hanging_nodes,CELL | FACE,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
983
        //cout << rank << ": After end " << std::endl;
984
		//reorder element's data to free up space
985
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
986
		m->ReorderEmpty(CELL|FACE|EDGE|NODE);
987
		EXIT_BLOCK();
988
989
		//return number of refined cells
		call_counter++;
Kirill Terekhov's avatar
Kirill Terekhov committed
990
		ret = m->Integrate(ret);
991
		REPORT_VAL("ret ",ret)
992
        EXIT_FUNC();
993
994
		return ret != 0;
	}
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006

    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;
    }
1007
	
1008

1009
1010
	bool AdaptiveMesh::Coarse(TagInteger & indicator)
	{
1011
		std::string file;
1012
		ENTER_FUNC();
Kirill Terekhov's avatar
Kirill Terekhov committed
1013
        //return false;
1014
1015
1016
1017
		static int call_counter = 0;
		//return number of coarsened cells
		int ret = 0;
		//initialize tree structure