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


#if defined(USE_PARALLEL_WRITE_TIME)
#define REPORT_MPI(x) {m->WriteTab(m->out_time) << "<MPI><![CDATA[" << #x << "]]></MPI>" << std::endl; x;}
#define REPORT_STR(x) {m->WriteTab(m->out_time) << "<TEXT><![CDATA[" << x << "]]></TEXT>" << std::endl;}
#define REPORT_VAL(str,x) {m->WriteTab(m->out_time) << "<VALUE name=\"" << str << "\"> <CONTENT><![CDATA[" << x << "]]></CONTENT> <CODE><![CDATA[" << #x << "]]></CODE></VALUE>" << std::endl;}
#define ENTER_FUNC() double all_time = Timer(); {m->WriteTab(m->out_time) << "<FUNCTION name=\"" << __FUNCTION__ << "\" id=\"func" << m->GetFuncID()++ << "\">" << std::endl; m->Enter();}
11
#define ENTER_BLOCK() { double btime = Timer(); m->WriteTab(m->out_time) << "<FUNCTION name=\"" << __FUNCTION__ << ":" << NameSlash(__FILE__) << ":" << __LINE__ << "\" id=\"func" << m->GetFuncID()++ << "\">" << std::endl; m->Enter();
12
13
14
15
16
17
18
#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
__INLINE std::string NameSlash(std::string input)
27
{
28
	for(unsigned l = static_cast<unsigned>(input.size()); l > 0; --l)
29
30
31
32
33
		if( input[l-1] == '/' || input[l-1] == '\\' )
			return std::string(input.c_str() + l);
	return input;
}

34
35
36
const bool check_orientation = false;
const bool check_convexity = false;

37

38
39
40
/// todo:
/// 1. coarsment
/// 2. strategy for faces/edges with faults
Kirill Terekhov's avatar
Kirill Terekhov committed
41
/// 3. geom model supportbn
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
/// 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));
	}
69
70
71
72
73
	void ReduceSum(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
	{
		(void) size;
		element.Real(tag) += *((const INMOST_DATA_REAL_TYPE *)data);
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
74
75
76
77
78
79
	
	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);
80
		std::vector<int> tmp(static_cast<size_t>(size)+static_cast<size_t>(odata.size()));
Kirill Terekhov's avatar
Kirill Terekhov committed
81
82
83
		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
84

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

98
        ss << " tag_owner (" << it.IntegerDF(m->OwnerTag()) << ")";
SilverLife's avatar
SilverLife committed
99
100
101

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

        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
115
            std::string type = "unknw";
SilverLife's avatar
SilverLife committed
116
117
118
119
            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
120
            ss << type << "-" << std::setw(2) << p->GlobalID() << " ";
SilverLife's avatar
SilverLife committed
121
122
            p++;
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
123
        ss << std::endl;
SilverLife's avatar
SilverLife committed
124
125
126
127
128
129

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

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

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

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

        if (face)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
184
            ss << "Faces: " << m->NumberOfFaces() <<  std::endl;
185
            for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) 
SilverLife's avatar
SilverLife committed
186
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
187
188
                ss << rank << ": " << std::setw(2) << it->LocalID() << " " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
189
190
191
192
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";

193
194
195

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

198
                ss << "  " << m->GetMarker(*it,m->NewMarker());
199
200
201
202

                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
203
                    ss << std::setw(2) << node->GlobalID() << " ";
204
                ss << ")";
Kirill Terekhov's avatar
Kirill Terekhov committed
205
                ss << std::endl;
SilverLife's avatar
SilverLife committed
206
207
208
209
210
            }
        }

        if (edge)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
211
            ss << "Edges: " << m->NumberOfEdges() <<  std::endl;
212
            for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it) 
SilverLife's avatar
SilverLife committed
213
214
215
216
217
            {
                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
218
                ss << std::endl;
SilverLife's avatar
SilverLife committed
219
220
221
222
223
            }
        }

        if (node)
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
224
            ss << "Nodes:" << std::endl;
225
            for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) 
SilverLife's avatar
SilverLife committed
226
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
227
228
                ss << rank << ": " << std::setw(2) << it->GlobalID() << " - " ;            
                ss << std::setw(6);
SilverLife's avatar
SilverLife committed
229
230
231
232
233
                if (it->GetStatus() == Element::Shared) ss << "shared";
                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
                else ss << "none";

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

                    ss << "(" << 
Kirill Terekhov's avatar
Kirill Terekhov committed
237
238
239
                        std::setw(3) << it->RealArray(m->CoordsTag())[0] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[1] << " " << 
                        std::setw(3) << it->RealArray(m->CoordsTag())[2] << ")";
240

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

Kirill Terekhov's avatar
Kirill Terekhov committed
246
247
        ss << "=========================================" << std::endl;
        os << ss.str() << std::endl;
SilverLife's avatar
SilverLife committed
248
    }
249
	*/
SilverLife's avatar
SilverLife committed
250

251
	/*
SilverLife's avatar
SilverLife committed
252
253
254
255
256
    void AdaptiveMesh::UpdateStatus()
    {

        for(ElementType mask = CELL; mask >= NODE; mask = PrevElementType(mask))
        {
257
            for(Mesh::iteratorElement it = m->BeginElement(mask); it != m->EndElement(); it++)
SilverLife's avatar
SilverLife committed
258
259
260
261
262
263
264
265
266
            {
                int stat = 0;
                if (it->GetStatus() == Element::Shared) stat = 1;
                else if (it->GetStatus() == Element::Ghost)  stat = 2;

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

285
	/*
286
287
288
    void AdaptiveMesh::SynchronizeSet(ElementSet set)
    {
    #ifdef USE_MPI
289
290
        int size = m->GetProcessorsNumber();
        int rank = m->GetProcessorRank();
291
292
        for (int i = 0; i < size; i++)
        {
293
294
            set.IntegerArray(m->SendtoTag()).push_back(i);
            m->ExchangeMarked();
295
296
297
        }
    #endif
    }
298
	*/
299

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

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

471
472
		//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
473
		
474
475
476
477
		//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();
478

479
		m->CheckSetLinks(__FILE__,__LINE__);
480
		CheckParentSet(__FILE__,__LINE__);
481
		/*
482
		ENTER_BLOCK();
483
		m->ResolveSets();
484
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
485
		m->ExchangeData(parent_set,CELL,0);
486
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
487
		m->ResolveSets();
488
		//m->CheckCentroids(__FILE__,__LINE__);
489
		//m->ExchangeData(parent_set,CELL,0);
490
		//m->CheckCentroids(__FILE__,__LINE__);
Kirill Terekhov's avatar
Kirill Terekhov committed
491
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
492
		//m->CheckCentroids(__FILE__,__LINE__);
493
		//CheckParentSet(Tag());
494
		EXIT_BLOCK();
495
		*/
496
497
498
499
		//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();
500

501
502
		//m->Save("before_refine_parent"+std::to_string(fi)+".pvtk");
		//std::cout << "Save before_refine_parent"+std::to_string(fi)+".pvtk" << std::endl;
503
504
505
		

		
506
		int schedule_counter = 1; //indicates order in which refinement will be scheduled
507
		Storage::integer scheduled = 1; //indicates that at least one element was scheduled on current sweep
508
		ENTER_BLOCK();
509
		//0. Extend indicator for edges and faces
Kirill Terekhov's avatar
Kirill Terekhov committed
510
		indicator = m->CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
511
512
		//1.Communicate indicator - it may be not synced
		m->ExchangeData(indicator,CELL,0);
513
514
		while(scheduled)
		{
515
			REPORT_VAL("scheduled",scheduled);
516
517
			//2.Propogate indicator down to the faces,edges
			//  select schedule for them
518
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
519
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
520
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
521
				Cell c = m->CellByLocalID(it);
522
523
524
525
526
527
528
529
530
531
				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
					}
				}
			}
532
			EXIT_BLOCK();
533
			//3.Communicate indicator on faces and edges
534
			m->ReduceData(indicator,FACE|EDGE,0,ReduceMax);
Kirill Terekhov's avatar
Kirill Terekhov committed
535
			m->ExchangeData(indicator,FACE|EDGE,0);
536
537
538
539
			//4.Check for each cell if there is
			//  any hanging node with adjacent in a need to refine,
			//  schedule for refinement earlier.
			scheduled = 0;
540
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
541
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
542
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
543
				Cell c = m->CellByLocalID(it);
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
				//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;
						}
					}
				}
			}
565
			EXIT_BLOCK();
566
567
568
569
			//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
570
			scheduled = m->Integrate(scheduled);
571
572
			if( scheduled ) schedule_counter++;
		}
573
		//m->ExchangeData(indicator,CELL | FACE | EDGE,0);
574
		EXIT_BLOCK();
575
576
		
		//m->ExchangeData(hanging_nodes,CELL|FACE,NONE);
Kirill Terekhov's avatar
Kirill Terekhov committed
577
		//m->Save("indicator.pmf");
578
579
580
581
582
583
		
		//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();
		
584
585
		m->ExchangeData(hanging_nodes,CELL | FACE,0);
		
586
587
		//if( !Element::CheckConnectivity(m) ) std::cout << __FILE__ << ":" << __LINE__ << " broken connectivity" << std::endl;
		//m->CheckCentroids(__FILE__,__LINE__);
588
		//6.Refine
589
590
591

		assert(Element::CheckConnectivity(m));
		CheckClosure(__FILE__,__LINE__);
592
		ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
593
		m->BeginModification();
594
595
		while(schedule_counter)
		{
596
597
			assert(Element::CheckConnectivity(m));
			CheckClosure(__FILE__,__LINE__);
598
599
			Storage::real xyz[3] = {0,0,0};
			//7.split all edges of the current schedule
600
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
601
			for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
602
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
603
				Edge e = m->EdgeByLocalID(it);
604
605
				if( !e.Hidden() && indicator[e] == schedule_counter )
				{
606
607
608
609
610
611
612
					//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 ") );
613
614
615
					//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
616
					for(Storage::integer d = 0; d < m->GetDimensions(); ++d)
617
618
619
						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
620
					Node n = m->CreateNode(xyz);
621
					//set increased level for new node
622
					level[n] = level[e.getBeg()] = level[e.getEnd()] = level[e]+1;
623
624
625
					//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);
626
					//CheckClosure(__FILE__,__LINE__);
627
					//split the edge by the middle node
Kirill Terekhov's avatar
Kirill Terekhov committed
628
					ElementArray<Edge> new_edges = Edge::SplitEdge(e,ElementArray<Node>(m,1,n.GetHandle()),0);
629
					for(ElementArray<Face>::size_type kt = 0; kt < edge_faces.size(); ++kt) assert(edge_faces[kt].Closure());
630
631
					//set increased level for new edges
					level[new_edges[0]] = level[new_edges[1]] = level[e]+1;
632
#if defined(USE_AUTODIFF) && defined(USE_SOLVER)
633
					if( model ) model->EdgeRefinement(e,new_edges);
634
#endif
635
636
637
638
					//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]]);
					//}
639
					//CheckClosure(__FILE__,__LINE__);
640
					//if( !Element::CheckConnectivity(m) ) std::cout << __FILE__ << ":" << __LINE__ << " broken connectivity" << std::endl;
641
642
				}
			}
643
			EXIT_BLOCK();
644
645
			assert(Element::CheckConnectivity(m));
			CheckClosure(__FILE__,__LINE__);
646
			//8.split all faces of the current schedule, using hanging nodes on edges
647
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
648
			for(Storage::integer it = 0; it < m->FaceLastLocalID(); ++it) if( m->isValidFace(it) )
649
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
650
				Face f = m->FaceByLocalID(it);
651
652
				if( !f.Hidden() && indicator[f] == schedule_counter )
				{
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
#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
689
690
691
692
693
694
695
696
697
698
699
700
701
					//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
702
					Node n = m->CreateNode(xyz);
703
					//set increased level for the new node
704
					level[n] = level[f]+1;
705
706
707
					//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
708
709
					ElementArray<Node> edge_nodes(m,2); //to create new edges
					ElementArray<Edge> hanging_edges(m,face_hanging_nodes.size());
710
711
712
713
					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
714
						hanging_edges[kt] = m->CreateEdge(edge_nodes).first;
715
716
717
718
719
720
721
722
						//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;
723
#if defined(USE_AUTODIFF) && defined(USE_SOLVER)
724
					if( model ) model->FaceRefinement(f,new_faces);
725
#endif
726
727
				}
			}
728
			EXIT_BLOCK();
729
730
			assert(Element::CheckConnectivity(m));
			CheckClosure(__FILE__,__LINE__);
731
			//this tag helps recreate internal face
Kirill Terekhov's avatar
Kirill Terekhov committed
732
			TagReferenceArray internal_face_edges = m->CreateTag("INTERNAL_FACE_EDGES",DATA_REFERENCE,NODE,NODE,4);
733
			//this marker helps detect edges of current cell only
Kirill Terekhov's avatar
Kirill Terekhov committed
734
			MarkerType mark_cell_edges = m->CreateMarker();
735
			//this marker helps detect nodes hanging on edges of unrefined cell
Kirill Terekhov's avatar
Kirill Terekhov committed
736
			MarkerType mark_hanging_nodes = m->CreateMarker();
737
			//9.split all cells of the current schedule
738
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
739
			for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
740
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
741
				Cell c = m->CellByLocalID(it);
742
743
744
745
746
747
748
749
750
751
752
				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
753
					Node n = m->CreateNode(xyz);
754
					//set increased level for the new node
755
					level[n] = level[c]+1;
756
757
					//retrive all edges of current face to mark them
					ElementArray<Edge> cell_edges = c.getEdges();
758
759
760
761
762
#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
763
764
765
					//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
766
767
					ElementArray<Node> edge_nodes(m,2);
					ElementArray<Edge> edges_to_faces(m,cell_hanging_nodes.size());
768
769
770
771
772
773
					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
774
						edges_to_faces[kt] = m->CreateEdge(edge_nodes).first;
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
						//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
810
					ElementArray<Face> internal_faces(m,edge_hanging_nodes.size());
811
812
813
814
815
816
817
818
819
820
					//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
821
						internal_faces[kt] = m->CreateFace(ElementArray<Edge>(m,face_edges.begin(),face_edges.end())).first;
822
823
824
825
826
						//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
827
828
829
830
831
832
833
					//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;
					//	
					//}
834
835
					//split the cell
					//retrive parent set
Kirill Terekhov's avatar
Kirill Terekhov committed
836
					ElementSet parent(m,parent_set[c]);
837
					//create set corresponding to old coarse cell
838
839
					//Storage::real cnt[3];
					//c.Centroid(cnt);
840
					std::stringstream set_name;
Kirill Terekhov's avatar
Kirill Terekhov committed
841
					//set_name << parent.GetName() << "_C" << c.GlobalID(); //rand may be unsafe
Kirill Terekhov's avatar
Kirill Terekhov committed
842
					if( parent == root )
843
						set_name << "AM_R" << set_id[c];
Kirill Terekhov's avatar
Kirill Terekhov committed
844
					else
845
						set_name << parent.GetName() << "C" << set_id[c];
846
					//set_name << base64_encode_((unsigned char *)cnt,3*sizeof(double)/sizeof(unsigned char));
Kirill Terekhov's avatar
Kirill Terekhov committed
847
					
Kirill Terekhov's avatar
Kirill Terekhov committed
848
849
850
851
852
853
					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;
854
						std::cout << rank << " Elements of " << check_set.GetName() << ": ";
Kirill Terekhov's avatar
Kirill Terekhov committed
855
						for(ElementSet::iterator it = check_set.Begin(); it != check_set.End(); ++it)
856
							std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << "," << indicator[*it] << " ";
Kirill Terekhov's avatar
Kirill Terekhov committed
857
						std::cout << std::endl;
858
859
860
861
862
863
864
865
866
867
868
						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
869
						exit(-1);
Kirill Terekhov's avatar
Kirill Terekhov committed
870
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
871
					
Kirill Terekhov's avatar
Kirill Terekhov committed
872
					ElementSet cell_set = m->CreateSetUnique(set_name.str()).first;
873
					//cell_set->SetExchange(ElementSet::SYNC_ELEMENTS_ALL);
874
					level[cell_set] = level[c]+1;
875
876
877
878
					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));
879
880
881
					//set up increased level for the new cells
					for(ElementArray<Cell>::size_type kt = 0; kt < new_cells.size(); ++kt)
					{
882
						set_id[new_cells[kt]] = (int)kt;
883
884
885
886
						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
887
888
889
890
891
892
893
894
895
896
					/*
					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() )
897
					parent.AddChild(cell_set);
898
#if defined(USE_AUTODIFF) && defined(USE_SOLVER)
899
					if( model ) model->CellRefinement(c,new_cells);
900
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
901
					//else assert(cell_set->GetParent() == parent);
902
903
904
905
					//increment number of refined cells
					ret++;
				}
			}
906
			EXIT_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
907
908
909
			m->ReleaseMarker(mark_hanging_nodes);
			m->ReleaseMarker(mark_cell_edges);
			m->DeleteTag(internal_face_edges);
910
911
			assert(Element::CheckConnectivity(m));
			CheckClosure(__FILE__,__LINE__);
912
913
914
			//10.jump to later schedule, and go to 7.
			schedule_counter--;
		}
915
		m->CheckSetLinks(__FILE__,__LINE__);
916
		//free created tag
Kirill Terekhov's avatar
Kirill Terekhov committed
917
		m->DeleteTag(indicator,FACE|EDGE);
SilverLife's avatar
SilverLife committed
918

919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
		if( check_orientation )
		{
			ENTER_BLOCK();
			int nfixed = 0, nfixednew = 0, nfixedbnd = 0, nghost = 0;
			for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) 
				if( !it->CheckNormalOrientation() )
				{
					//it->FixNormalOrientation();
					nfixed++;
					if( it->New() )
						nfixednew++;
					if( it->Boundary() )
						nfixedbnd++;
					if( it->GetStatus() == Element::Ghost )
						nghost++;
					std::cout << "rank " << m->GetProcessorRank() << " face " << it->GetHandle() << std::endl;
				}
				//std::cout << "Face " << it->LocalID() << " oriented incorrectly " << std::endl;
			if( nfixed ) 
			{
				std::cout << __FILE__ << ":" << __LINE__ << " rank " << rank << " bad " << nfixed << " (new " << nfixednew << " bnd " << nfixedbnd << " ghost " << nghost << ") faces " << std::endl;
				REPORT_STR(rank << " bad " << nfixed << " faces");
			}
			EXIT_BLOCK();
		}
		
		if( check_convexity )
		{
			int nbad = 0;
			for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
				if( !it->CheckConvexity() ) nbad++;
			if( nbad ) std::cout << __FILE__ << ":" << __LINE__ << " rank " << rank << " nonconvex cells: " << nbad << std::endl;
		}

SilverLife's avatar
SilverLife committed
953

954
		m->CheckSetLinks(__FILE__,__LINE__);
955
		//11. Restore parallel connectivity, global ids
Kirill Terekhov's avatar
Kirill Terekhov committed
956
		m->ResolveModification();
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990

		if( check_orientation )
		{
			ENTER_BLOCK();
			int nfixed = 0, nfixednew = 0, nfixedbnd = 0, nghost = 0;
			for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) 
				if( !it->CheckNormalOrientation() )
				{
					//it->FixNormalOrientation();
					nfixed++;
					if( it->New() )
						nfixednew++;
					if( it->Boundary() )
						nfixedbnd++;
					if( it->GetStatus() == Element::Ghost )
						nghost++;
					std::cout << "rank " << m->GetProcessorRank() << " face " << it->GetHandle() << std::endl;
				}
				//std::cout << "Face " << it->LocalID() << " oriented incorrectly " << std::endl;
			if( nfixed ) 
			{
				std::cout << __FILE__ << ":" << __LINE__ << " rank " << rank << " bad " << nfixed << " (new " << nfixednew << " bnd " << nfixedbnd << " ghost " << nghost << ") faces " << std::endl;
				REPORT_STR(rank << " bad " << nfixed << " faces");
			}
			EXIT_BLOCK();
		}
		
		if( check_convexity )
		{
			int nbad = 0;
			for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
				if( !it->CheckConvexity() ) nbad++;
			if( nbad ) std::cout << __FILE__ << ":" << __LINE__ << " rank " << rank << " nonconvex cells: " << nbad << std::endl;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
991
		//m->SynchronizeMarker(m->NewMarker(),CELL|FACE|EDGE|NODE,SYNC_BIT_OR);
992
        //ExchangeGhost(3,NODE); // Construct Ghost cells in 2 layers connected via nodes
993
		//12. Let the user update their data
994
		//todo: call back function for New( cells
995
#if defined(USE_AUTODIFF) && defined(USE_SOLVER)
Kirill Terekhov's avatar
Kirill Terekhov committed
996
		if( model ) model->Adaptation(*m);
Kirill Terekhov's avatar