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

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

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

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

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

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

955
		m->CheckSetLinks(__FILE__,__LINE__);
956
		//11. Restore parallel connectivity, global ids
Kirill Terekhov's avatar
Kirill Terekhov committed
957
		m->ResolveModification();
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
991

		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
992
		//m->SynchronizeMarker(m->NewMarker(),CELL|FACE|EDGE|NODE,SYNC_BIT_OR);
993
        //ExchangeGhost(3,NODE); // Construct Ghost cells in 2 layers connected via nodes
994
		//12. Let the user update their data
995
		//todo: call back function for New( cells
996
#if defined(USE_AUTODIFF) && defined(USE_SOLVER)
Kirill Terekhov's avatar