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

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

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

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

						ElementArray<Cell> new_cells = Cell::SplitCell(c,internal_faces,0);
						splits++;
						t2 = Timer(), tsplit += t2 - t1, t1 = t2;
						
						std::sort(new_cells.begin(),new_cells.end(),Mesh::CentroidComparator(m));
						
						t2 = Timer(), tsort += t2 - t1, t1 = t2;
						//set up increased level for the new cells
						for(ElementArray<Cell>::size_type kt = 0; kt < new_cells.size(); ++kt)
954
						{
955
956
957
958
959
960
961
962
963
964
965
							set_id[new_cells[kt]] = (int)kt;
							level[new_cells[kt]] = level[c]+1;
							cell_set.PutElement(new_cells[kt]);
							parent_set[new_cells[kt]] = cell_set.GetHandle();
						}
						/*
						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] << " ";
966
967
							std::cout << std::endl;
						}
968
969
970
971
						*/
						//if( !cell_set->HaveParent() )
						parent.AddChild(cell_set);
						t2 = Timer(), tpconn += t2 - t1, t1 = t2;
972
#if defined(USE_AUTODIFF) && defined(USE_SOLVER)
973
						if( model ) model->CellRefinement(c,new_cells);
974
#endif
975
976
977
978
979
						t2 = Timer(), tdata += t2 - t1, t1 = t2;
						//else assert(cell_set->GetParent() == parent);
						//increment number of refined cells
						ret++;
					}
980
				}
981
982
983
984
985
986
987
988
989
990
991
992
993
994
				REPORT_VAL("adjacencies", tadj);
				REPORT_VAL("create (new nodes)", tcreate);
				REPORT_VAL("hanging1 (new edges)", thanging1);
				REPORT_VAL("hanging2 (new faces)", thanging2);
				REPORT_VAL("create set", tset);
				REPORT_VAL("connect set", tpconn);
				REPORT_VAL("data", tdata);
				REPORT_VAL("split", tsplit);
				REPORT_VAL("sort cells", tsort);
				REPORT_VAL("new nodes", new_nodes);
				REPORT_VAL("new edges", new_edges);
				REPORT_VAL("new faces", new_faces);
				REPORT_VAL("new sets", new_sets);
				REPORT_VAL("splits", splits);
995
			}
996
			EXIT_BLOCK();
997
			ENTER_BLOCK();
Kirill Terekhov's avatar
Kirill Terekhov committed
998
999
1000
			m->ReleaseMarker(mark_hanging_nodes);
			m->ReleaseMarker(mark_cell_edges);
			m->DeleteTag(internal_face_edges);
1001
1002
			EXIT_BLOCK();
			ENTER_BLOCK();
1003
1004
			assert(Element::CheckConnectivity(m));
			CheckClosure(__FILE__,__LINE__);
1005
			EXIT_BLOCK();
1006
			//10.jump to later schedule, and go to 7.
1007
			REPORT_VAL("schedule counter",schedule_counter);
1008
1009
			schedule_counter--;
		}
1010
		m->CheckSetLinks(__FILE__,__LINE__);
1011
		//free created tag
Kirill Terekhov's avatar
Kirill Terekhov committed
1012
		m->DeleteTag(indicator,FACE|EDGE);
SilverLife's avatar
SilverLife committed
1013

1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
		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 )
		{
1042
			ENTER_BLOCK();
1043
1044
1045
1046
			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;
1047
			EXIT_BLOCK();
1048
		}
1049
		ENTER_BLOCK();
1050
		m->CheckSetLinks(__FILE__,__LINE__);
1051
1052
1053
1054
		EXIT_BLOCK();
		ENTER_BLOCK();
		m->Barrier();
		EXIT_BLOCK();
1055
		//11. Restore parallel connectivity, global ids
Kirill Terekhov's avatar
Kirill Terekhov committed
1056
		m->ResolveModification();
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085