mesh_ecl_file.cpp 188 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
#ifdef _MSC_VER //kill some warnings
#define _CRT_SECURE_NO_WARNINGS
Kirill Terekhov's avatar
Kirill Terekhov committed
3
#define _SCL_SECURE_NO_WARNINGS
Kirill Terekhov's avatar
Kirill Terekhov committed
4
5
6
#endif

#include "inmost.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
7
#include "../Mesh/incident_matrix.hpp"
Kirill Terekhov's avatar
Kirill Terekhov committed
8
#include <string>
Kirill Terekhov's avatar
Kirill Terekhov committed
9
#include <ctime>
Kirill Terekhov's avatar
Kirill Terekhov committed
10
11
#if defined(USE_MESH)

Kirill Terekhov's avatar
Kirill Terekhov committed
12
13
14
15
16
17
18
19
20
21
// coords/zcorn algorithm
// 1) put all block nodes into pillars, sort each pillar nodes by actual depth
// 2) create edges down along pillars and mark them according to blocks
// 3) consider all pairs of pillars in Oxz and Oyz planes and create uncut block edges, mark them in block numbers
// 4) use line-sweep algorithm to intersect and cut all the edges between each pair of pillars, mark them in block numbers
// 5) mark all the nodes with block numbers that they belong by considering union of block numbers on adjacent edges
// 5) use incident_matrix algorithm to create all the faces between pair of pillars from pillar edges and inter-pillar edges
// 6) from intersection of block numbers on nodes figure out which blocks the face belongs
// 7) add top and bottom interface

Kirill Terekhov's avatar
Kirill Terekhov committed
22
//todo
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
23
24
// 1. (ok) remove edges on pillars that do not get any associated block
// 2. (ok) do not create cells with actnum = 0
Kirill Terekhov's avatar
Kirill Terekhov committed
25
// 2. fix fast intersect() algorithm, check against intersect_naive()
26
27
28
29
// 3. (ok) CheckEdgeOrder reports bad order of edges for certain interfaces - investigate
// 4. (ok) when encounter cells with one face with actnum = 1 should disconnect adjacent cells, see TODO4: in text
// 4.1 (ok) with TRANM
// 4.2 (ok) with Fracture algorithm
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
30
// 5. (ok) populate keywords data, poro, perm, actnum, satnum, etc...
31
32
33
// 5.1 (ok) add saturation, pressure
// 5.2 (ok) eqlnum, rocknum, pvtnum
// 5.3 add xmf, ymf, zmf, read number of components
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
34
// 6. local grid refinement
35
36
// 6.1 radial local grid refinement
// 6.2 nested local grid refinement
37
// 7. (ok) omp-parallel loading
38
// 8. (ok) mpi-parallel loading
39
// 9. (ok) do not convert arrays xyz and zcorn
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
40
//10. read wells into sets
Kirill Terekhov's avatar
Kirill Terekhov committed
41
//10.1 (ok) read compdat
Kirill Terekhov's avatar
Kirill Terekhov committed
42
43
44
45
//10.2 (ok) read properties
//10.3 RESV, TODO10.3
//10.4 schedule for wells
//10.5 track previous record on schedule to optmize number of records
46
47
//11. optimize:
//11.1 detect that pillars have no faults, skip intersection in this case
48
49
50
51
//11.2 (ok) skip incident_matrix algorithm when no inner edges found
//11.3 drop out edges from intersection that are not expected to intersect
//11.4 if intersections array is empty consider faster algorithm
//11.5 (slower) tree along z direction in transformed coordinates in intersect_naive
Kirill Terekhov's avatar
Kirill Terekhov committed
52
53
//12. (cancel, no longer needed with new method for intersection point position) High-order elements: edge, face for ECL_CURVILINEAR option, see TODO12, need high order elements support and algorithms in geometry.cpp
//13. (fail, projection onto triangles is very bad) When curvature is high, split edges along triangle edge before intersection, see TODO13
54
55
56
57
//14. Process BOX/ENDBOX keywords
//15. (ok) verbosity
//16. Make cell centers inside blocks TODO16
//17. (ok) Erase unused tags
Kirill Terekhov's avatar
Kirill Terekhov committed
58

59
#define ECLSTRCMP(x,y) strncmp(STR8(x).c_str(),STR8(y).c_str(),8)
Kirill Terekhov's avatar
Kirill Terekhov committed
60

Kirill Terekhov's avatar
Kirill Terekhov committed
61
62
63
64
65
66

//controls when to consider two nodes on pillar to be the same
#define ECL_PILLAR_DEPTH_EPS 1.0e-8
#define ECL_POINT_EPS 1.0e-8
#define ECL_INTERSECT_EPS 1.0e-12

Kirill Terekhov's avatar
Kirill Terekhov committed
67
68
69
//this controls what is recorded into tag CTRL of well sets
#define ECL_WCTRL_RATE 0
#define ECL_WCTRL_BHP 1
70
#define ECL_WCTRL_RESV 2
Kirill Terekhov's avatar
Kirill Terekhov committed
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89

//this controls what is recorded into tag PHASE of well sets
#define ECL_WPHASE_WATER 0
#define ECL_WPHASE_OIL 1
#define ECL_WPHASE_GAS 3
#define ECL_WPHASE_MULTI 4 //LRAT

//this controls what is recorded into tag TYPE of well sets
#define ECL_WTYPE_INJE 0
#define ECL_WTYPE_PROD 1


//this controls what tag is changed by schedule record
#define ECL_WTAG_WI 0
#define ECL_WTAG_CVAL 1
#define ECL_WTAG_STATE 3
#define ECL_WTAG_CTRL 4
#define ECL_WTAG_TYPE 5
#define ECL_WTAG_PHASE 6
90
#define ECL_WTAG_BHP 7
Kirill Terekhov's avatar
Kirill Terekhov committed
91

92
//eclipse states
Kirill Terekhov's avatar
Kirill Terekhov committed
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#define ECL_NEVER -1
#define ECL_NONE 0
#define ECL_SKIP_SECTION 1
#define ECL_INCLUDE 2
#define ECL_DIMENS 3
#define ECL_DX 4
#define ECL_DY 5
#define ECL_DZ 6
#define ECL_TOPS 7
#define ECL_PERMX 8
#define ECL_PERMY 9
#define ECL_PERMZ 10
#define ECL_PORO 11
#define ECL_MAPAXIS 12
#define ECL_INRAD 13
#define ECL_COORDS 14
#define ECL_ZCORN 15
110
#define ECL_ACTNUM 16
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
111
#define ECL_SATNUM 17
Kirill Terekhov's avatar
Kirill Terekhov committed
112
#define ECL_THCONR 18
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#define ECL_MULTIPLY 19
#define ECL_MULTIPLY_MUL 20
#define ECL_MULTIPLY_BLK 21
#define ECL_EDITNNC 22
#define ECL_EDITNNC_BLK 23
#define ECL_EDITNNC_MUL 24
#define ECL_NTG 25
#define ECL_FAULTS 26
#define ECL_FAULTS_BLK 27
#define ECL_FAULTS_DIR 28
#define ECL_MULTFLT 29
#define ECL_MULTFLT_MUL 30
#define ECL_PRESSURE 31
#define ECL_SWAT 32
#define ECL_SOIL 33
#define ECL_SGAS 34
#define ECL_PBUB 35
#define ECL_EQLNUM 36
#define ECL_ROCKNUM 37
#define ECL_PVTNUM 38
Kirill Terekhov's avatar
Kirill Terekhov committed
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
#define ECL_COMPDAT 39
#define ECL_COMPDAT_BLK 40
#define ECL_COMPDAT_OPEN 41
#define ECL_COMPDAT_SATNUM 42
#define ECL_COMPDAT_TRANS 43
#define ECL_COMPDAT_BORE 44
#define ECL_COMPDAT_PERM 45
#define ECL_COMPDAT_SKIN 46
#define ECL_COMPDAT_DFAC 47
#define ECL_COMPDAT_DIR 48
#define ECL_COMPDAT_RAD 49
#define ECL_DATES 50
#define ECL_DATES_MON 51
#define ECL_DATES_YEAR 52
#define ECL_DATES_HRS 53
#define ECL_WELSPECS 54
#define ECL_WELSPECS_GROUP 55
#define ECL_WELSPECS_I 56
#define ECL_WELSPECS_J 57
#define ECL_WELSPECS_Z 58
#define ECL_WELSPECS_PHASE 59
#define ECL_TSTEP 60
#define ECL_WCONPRODINJE 61
#define ECL_WCONPRODINJE_TYPE 62
#define ECL_WCONPRODINJE_OPEN 63
#define ECL_WCONPRODINJE_CTRL 64
#define ECL_WCONPRODINJE_RATES 65
#define ECL_WCONPRODINJE_BHP 66
161

Kirill Terekhov's avatar
Kirill Terekhov committed
162
163
164
165
166
167
168
169
170
171
172
173

#define ECL_GTYPE_NONE 0
#define ECL_GTYPE_TOPS 1
#define ECL_GTYPE_ZCORN 2
#define ECL_GTYPE_RADIAL 3
#define ECL_GTYPE_CARTESIAN 4


#define ECL_VAR_NONE 0
#define ECL_VAR_REAL 1
#define ECL_VAR_INT 2

174
#define ECL_IJK_DATA(i,j,k) (i + ((j)+(k)*dims[1])*dims[0])
175
176
#define ECL_IJK_COORDS(i,j,l,coord) ((((i) + (j)*(dims[0]+1))*2+(l))*3+(coord))
#define ECL_IJK_ZCORN(i,j,k,l) ((k)*dims[1]*dims[0]*8 + (1-(l)/4)*dims[1]*dims[0]*4 + (j)*dims[0]*4 + ((l)/2%2)*dims[0]*2 + (i)*2 + (l)%2)
177
178
179
180
181

//line sweep events
#define SEG_START 0
#define SEG_END 1

Kirill Terekhov's avatar
Kirill Terekhov committed
182
183
184
#define HAVE_PERM_X 1
#define HAVE_PERM_Y 2
#define HAVE_PERM_Z 4
Kirill Terekhov's avatar
Kirill Terekhov committed
185

186

187

Kirill Terekhov's avatar
Kirill Terekhov committed
188
189
namespace INMOST
{
190

Kirill Terekhov's avatar
Kirill Terekhov committed
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269

	static void GetDXYZ(Element e, double & dx, double & dy, double & dz)
	{
		Storage::real maxmin[6];
		maxmin[0] = -1e20;
		maxmin[1] = 1e20;
		maxmin[2] = -1e20;
		maxmin[3] = 1e20;
		maxmin[4] = -1e20;
		maxmin[5] = 1e20;
		ElementArray<Node> nodes = e->getNodes();
		for (ElementArray<Node>::iterator it = nodes.begin(); it != nodes.end(); it++)
		{
			Storage::real_array c = it->Coords();
			for (int i = 0; i < (int)c.size(); i++)
			{
				if (maxmin[2 * i + 0] < c[i]) maxmin[2 * i + 0] = c[i]; //max
				if (maxmin[2 * i + 1] > c[i]) maxmin[2 * i + 1] = c[i]; //min
			}
		}
		dx = maxmin[0] - maxmin[1];
		dy = maxmin[2] - maxmin[3];
		dz = maxmin[4] - maxmin[5];
	}

	bool ReadName(const char * p, const char * pend, char * rec, int * nchars)
	{
		int k = 0;
		int r = 0;
		int s;
		while (p + k != pend && isspace(p[k])) ++k;
		s = k;
		if (p[k] == '\'' || p[k] == '"')
		{
			char quote = p[k++];
			while (p + k != pend && p[k] != quote)
				rec[r++] = p[k++];
			*nchars = k+1; //skip quote
		}
		else
		{
			while (p + k != pend && !isspace(p[k]))
				rec[r++] = p[k++];
			*nchars = k;
		}
		rec[r] = '\0';
		return s != k;
	}

	int MonthInt(std::string mo)
	{
		if( mo == "JAN" )
			return 1;
		else if( mo == "FEB" )
			return 2;
		else if( mo == "MAR" )
			return 3;
		else if( mo == "APR" )
			return 4;
		else if( mo == "MAY" )
			return 5;
		else if( mo == "JUN" )
			return 6;
		else if( mo == "JLY" )
			return 7;
		else if( mo == "JUL" )
			return 7;
		else if( mo == "AUG" )
			return 8;
		else if( mo == "SEP" )
			return 9;
		else if( mo == "OCT" )
			return 10;
		else if( mo == "NOV" )
			return 11;
		else if( mo == "DEC" )
			return 12;
		return -1;
	}
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
	void Faults(Mesh * m, MarkerType frac_markers)
	{
		//m->BeginModification();
		const Storage::real mult = 0.95;
		Tag indexes = m->CreateTag("FRAC_TEMP_INDEXES", DATA_INTEGER, CELL, NONE, 2);
		//associate a node or a pair of nodes to each non-fracture edge that ends with a fracture node
		Tag cell2node = m->CreateTag("CELL2NODE", DATA_REFERENCE, CELL, NONE);
		//new faces to reconnect
		Tag connface = m->CreateTag("CONNFACE", DATA_REFERENCE, CELL, NONE);

		//Unmark any boundary faces
		for (Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it)
		{
			if (!it->FrontCell().isValid() && it->GetMarker(frac_markers))
				it->RemMarker(frac_markers);
		}

		//For each node create several new ones that will be used in open fracture
		for (Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it)
		{
			ElementArray<Face> n_faces = it->getFaces(); //retrive all faces joining at the node
			INMOST_DATA_ENUM_TYPE num_frac_faces = 0;
			for (int k = 0; k < n_faces.size(); ++k) if (n_faces[k].GetMarker(frac_markers)) num_frac_faces++;
			if (num_frac_faces)
			{
				ElementArray<Cell> n_cells = it->getCells(); //get all cells of the node
				for (int k = 0; k < n_cells.size(); ++k) //assign local enumeration to the cells
					n_cells[k].IntegerDF(indexes) = k;
				for (int k = 0; k < n_faces.size(); ++k) //stich together node's numbers corresponding to cells if no fracture separates them
				{
					if (!n_faces[k].GetMarker(frac_markers) && n_faces[k].FrontCell().isValid())
					{
						int bi = n_faces[k].BackCell()->IntegerDF(indexes);
						int fi = n_faces[k].FrontCell()->IntegerDF(indexes);
						if (bi != fi)
						{
							for (int q = 0; q < n_cells.size(); q++) if (n_cells[q].IntegerDF(indexes) == fi)
								n_cells[q].IntegerDF(indexes) = bi;
						}
					}
				}
				dynarray<int, 64> nums(n_cells.size()); //store all numbers
				for (int k = 0; k < n_cells.size(); ++k) nums[k] = n_cells[k].IntegerDF(indexes);
				std::sort(nums.begin(), nums.end());
				nums.resize(std::unique(nums.begin(), nums.end()) - nums.begin());
				if (nums.size() > 1) //at least two distinctive nodes
				{
					for (int k = 0; k < nums.size(); k++)
					{
						Node image;
						Storage::real xyz[3] = { 0, 0, 0 }, cntc[3] = { 0, 0, 0 };
						int n_node_cells = 0;
						for (int q = 0; q < n_cells.size(); q++)
						{
							if (n_cells[q].IntegerDF(indexes) == nums[k])
							{
								n_cells[q].Centroid(cntc);
								xyz[0] += cntc[0];
								xyz[1] += cntc[1];
								xyz[2] += cntc[2];
								n_node_cells++;
							}
						}
						xyz[0] /= static_cast<Storage::real>(n_node_cells);
						xyz[1] /= static_cast<Storage::real>(n_node_cells);
						xyz[2] /= static_cast<Storage::real>(n_node_cells);
						Storage::real_array coords = it->Coords();
						for (Storage::real_array::size_type q = 0; q < coords.size(); ++q)
							xyz[q] = coords[q] * mult + (1 - mult)*xyz[q];
						image = m->CreateNode(xyz);

						for (int q = 0; q < n_cells.size(); q++)
						{
							if (n_cells[q].IntegerDF(indexes) == nums[k])
							{
								Storage::reference_array fnodes = n_cells[q].ReferenceArray(cell2node);
								fnodes.push_back(*it); //TODO
								fnodes.push_back(image);
							}
						}
					}
					it->SetMarker(frac_markers); //mark node as fracture node
					//std::cout << "Number of nodes: " << m->NumberOfNodes() << std::endl;
				}

			}
		}
		m->DeleteTag(indexes);


		//Mark edges that will be incorporated into fracture cells
		for (Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
		{
			if (it->nbAdjElements(FACE, frac_markers) > 2) it->SetMarker(frac_markers);
		}
		//create new matrix-matrix faces, adjacent to fractures
		//for all non-fracture faces that have any fracture edge create new face
		for (Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if (!it->GetMarker(frac_markers) && it->nbAdjElements(NODE, frac_markers))
		{
			ElementArray<Edge> edges = it->getEdges(); //go through edges of current face
			ElementArray<Edge> newedges(m);
			Storage::reference_array images = it->BackCell()->ReferenceArray(cell2node);
			//mark me and all my edges that they may contain images
			for (ElementArray<Edge>::iterator jt = edges.begin(); jt != edges.end(); ++jt)
			{
				if (jt->getBeg()->GetMarker(frac_markers) || jt->getEnd()->GetMarker(frac_markers))
				{
					ElementArray<Node> enodes(m, 2);
					if (jt->getBeg()->GetMarker(frac_markers)) //search image of node between my elements
					{
						for (int k = 0; k < static_cast<int>(images.size()) && !enodes[0].isValid(); k += 2)
						{
							if (images[k] == jt->getBeg()) // find my element between images
								enodes[0] = images[k + 1]->getAsNode();
						}
						assert(enodes[0] != InvalidElement());
					}
					else enodes[0] = jt->getBeg();
					if (jt->getEnd()->GetMarker(frac_markers))
					{
						for (int k = 0; k < static_cast<int>(images.size()) && !enodes[1].isValid(); k += 2)
						{
							if (images[k] == jt->getEnd()) // find my element between images
								enodes[1] = images[k + 1]->getAsNode();
						}
						assert(enodes[1] != InvalidElement());
					}
					else enodes[1] = jt->getEnd();
					if (enodes[0] != enodes[1]) newedges.push_back(m->CreateEdge(enodes).first);
				}
				else newedges.push_back(*jt);
			}
			std::pair<Face, bool> f = m->CreateFace(newedges);
			//indicate that this new face should be connected to respecting cells
			if (it->BackCell().isValid()) it->BackCell()->ReferenceArray(connface).push_back(f.first);
			if (it->FrontCell().isValid()) it->FrontCell()->ReferenceArray(connface).push_back(f.first);
		}

		//now create fracture-fracture control volumes, add fracture-matrix faces to matrix cells
		for (Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if (it->GetMarker(frac_markers) && it->nbAdjElements(NODE, frac_markers))
		{
			//I need all elements adjacent to neighbouring cells of current face
			//and ending at any element of current face
			//ElementArray<Face> fracfaces(m);
			ElementArray<Node> nodes = it->getNodes(), nodesb(m, nodes.size()), nodesf(m, nodes.size());
			ElementArray<Edge> edges = it->getEdges(), edgesb(m, nodes.size()), edgesf(m, nodes.size());
			INMOST_DATA_ENUM_TYPE N = (INMOST_DATA_ENUM_TYPE)nodes.size();
			Storage::reference_array images;
			{
				Cell c = it->BackCell();
				images = c->ReferenceArray(cell2node);
				//among images of fracture nodes select those highlighted by marker
				for (int q = 0; q < nodes.size(); ++q) if (nodes[q].GetMarker(frac_markers))
				{
					Node n;
					for (int k = 0; k < (int)images.size() && !n.isValid(); k += 2)
					if (images[k] == nodes[q]) n = images[k + 1]->getAsNode();
					assert(n.isValid());
					nodesb[q] = n;
				}
				else nodesb[q] = nodes[q];
				//create edges
				for (int k = 0; k < edges.size(); k++)
				{
					if (nodes[k].GetMarker(frac_markers) || nodes[(k + 1) % N].GetMarker(frac_markers))
					{
						ElementArray<Node> enodes(m, 2);
						enodes[0] = nodesb[k];
						enodes[1] = nodesb[(k + 1) % N];
						edgesb[k] = m->CreateEdge(enodes).first;
					}
					else edgesb.data()[k] = edges.data()[k];
				}
				//This is matrix-fracture face
				std::pair<Face, bool> facesb = m->CreateFace(edgesb);
				//add faces to indicate reconnection
				it->BackCell()->ReferenceArray(connface).push_back(facesb.first);
			}

			if (it->FrontCell().isValid()) //there is another cell
			{
				Cell c = it->FrontCell();
				images = c->ReferenceArray(cell2node);
				//among images of fracture nodes select those highlited by marker
				for (int q = 0; q < nodes.size(); ++q) if (nodes[q].GetMarker(frac_markers))
				{
					Node n;
					for (int k = 0; k < (int)images.size() && !n.isValid(); k += 2)
					if (images[k] == nodes[q]) n = images[k + 1]->getAsNode();
					assert(n.isValid());
					nodesf[q] = n;
				}
				else nodesf[q] = nodes[q];
				//create edges
				for (int k = 0; k < edges.size(); k++)
				{
					if (nodes[k].GetMarker(frac_markers) || nodes[(k + 1) % N].GetMarker(frac_markers))
					{
						ElementArray<Node> enodes(m, 2);
						enodes[0] = nodesf[k];
						enodes[1] = nodesf[(k + 1) % N];
						edgesf[k] = m->CreateEdge(enodes).first;
					}
					else edgesf[k] = edges[k];
				}
				//This is matrix-fracture face
				std::pair<Face, bool> facesf = m->CreateFace(edgesf);
				//add faces to indicate reconnection
				it->FrontCell()->ReferenceArray(connface).push_back(facesf.first);
			}
		}
		//now reconnect matrix cells to new faces
		for (Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) if (!it->GetMarker(frac_markers) && it->nbAdjElements(NODE, frac_markers))
		{
			ElementArray<Face> change(m);
			ElementArray<Face> faces = it->getFaces();
			for (ElementArray<Face>::iterator jt = faces.begin(); jt != faces.end(); jt++)
			{
				if (jt->nbAdjElements(NODE, frac_markers))
					change.push_back(*jt);
			}
			if (!change.empty())
			{
				Storage::reference_array newfaces = it->ReferenceArray(connface);
				assert(change.size() == newfaces.size());
				it->Disconnect(change.data(), (Storage::enumerator)change.size());
				it->Connect(newfaces.data(), (Storage::enumerator)newfaces.size());
			}
		}

		m->DeleteTag(connface);
		m->DeleteTag(cell2node);
		for (Mesh::iteratorElement it = m->BeginElement(FACE | EDGE | NODE); it != m->EndElement(); ++it)
		{
			it->RemMarker(frac_markers);
		}
		//delete remenants
		for (ElementType etypei = CELL; etypei > NODE; etypei = etypei >> 1)
		{
			ElementType etype = etypei >> 1;
			for (Mesh::iteratorElement it = m->BeginElement(etype); it != m->EndElement(); ++it)
			{
				if (!it->nbAdjElements(etypei))
					it->Delete();
			}
		}
		//m->ApplyModification();
		//m->EndModification();
		//m->Save("fracture_test.vtk");
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
	//replace n* witn n stars
	static void UnStar(const char * input, char * output)
	{
		size_t q = 1, j;
		output[0] = input[0]; //we skip first char
		if (input[0] != '\0')
		{
			for (size_t k = 1; k < strlen(input) - 1; ++k)
			{
				//current character is star, next char is space or end of line and prev char is not space
				if (input[k] == '*' && (isspace(input[k + 1]) || input[k + 1] == '/') && !isspace(input[k - 1]))
				{
					j = k;
					//roll back to previous space
					while (!isspace(input[k]) && k > 0)
					{
						k--;
						q--;
					}
					//advance from space
					if (isspace(input[k]))
					{
						k++;
						q++;
					}
					int count = atoi(input+k);
					for (int l = 0; l < count; ++l)
					{
						output[q++] = '*';
						output[q++] = ' ';
					}
					k = j;
				}
				else output[q++] = input[k];
			}
			output[q++] = input[strlen(input) - 1]; //we skip last char
			output[q++] = '\0'; // end of line
		}
	}

561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
	static std::string STR8(const char * input)
	{
		std::string ret(8, ' ');
		int kend = (int)strlen(input);
		if (kend > 8) kend = 8;
		for (int k = 0; k < kend; k++) ret[k] = input[k];
		return ret;
	}

	static std::string UnQoute(std::string input)
	{
		if (input[0] == '\'' && input[input.size() - 1] == '\'')
			return input.substr(1, input.size() - 2);
		else if (input[0] == '"' && input[input.size() - 1] == '"')
			return input.substr(1, input.size() - 2);
		return input;
	}
	static std::string ToUpper(std::string input)
	{
		for (size_t k = 0; k < input.size(); ++k)
			input[k] = toupper(input[k]);
		return input;
	}

	static int GetDir(std::string dir)
	{
		if (dir == "X-" || dir == "I-")
			return 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
589
		else if (dir == "X" || dir == "I" || dir == "X+" || dir == "I+" || dir == "X " || dir == "I ")
590
591
592
			return 1;
		else if (dir == "Y-" || dir == "J-")
			return 2;
Kirill Terekhov's avatar
Kirill Terekhov committed
593
		else if (dir == "Y" || dir == "J" || dir == "Y+" || dir == "J+" || dir == "Y " || dir == "J ")
594
595
596
			return 3;
		else if (dir == "Z-" || dir == "K-")
			return 4;
Kirill Terekhov's avatar
Kirill Terekhov committed
597
		else if (dir == "Z" || dir == "K" || dir == "Z+" || dir == "K+" || dir == "Z " || dir == "K ")
598
599
600
			return 5;
		return -1;
	}
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
601
602
603
	static std::string GetFolder(std::string file)
	{
		size_t found = file.find_last_of("/\\");
604
		if (found == std::string::npos)
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
605
			return "";
606
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
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
		else return file.substr(0, found);
	}

	static ElementArray<Edge> OrderEdges(const ElementArray<Edge> & input, MarkerType mrk)
	{
		Mesh * m = const_cast<Mesh *>(input.GetMeshLink());
		ElementArray<Edge> ret(m);
		ret.reserve(input.size());
		HandleType last;
		ret.push_back(input[0]);
		HandleType add = InvalidHandle();
		const Element::adj_type & lc0 = m->LowConn(input.at(0));
		for (int l = 1; l < (int)input.size(); ++l)
		{
			const Element::adj_type & lcl = m->LowConn(input.at(l));
			if (lcl[0] == lc0[0] || lcl[0] == lc0[1])
			{
				add = input.at(l);
				last = lcl[1];
				break;
			}
			else if (lcl[1] == lc0[0] || lcl[1] == lc0[1])
			{
				add = input.at(l);
				last = lcl[0];
				break;
			}
		}
		if (add == InvalidHandle())
			std::cout << __FILE__ << ":" << __LINE__ << " bad edges, number of edges " << input.size() << std::endl;
		else
		{
			ret.push_back(add);
			m->SetPrivateMarker(add, mrk);
		}

		while (ret.size() != input.size())
		{
			add = InvalidHandle();
			for (int l = 1; l < (int)input.size(); ++l)
			{
				if (!m->GetPrivateMarker(input.at(l), mrk))
				{
					const Element::adj_type & lcl = m->LowConn(input.at(l));
					if (lcl[0] == last)
					{
						add = input.at(l);
						last = lcl[1];
						break;
					}
					else if (lcl[1] == last)
					{
						add = input.at(l);
						last = lcl[0];
						break;
					}
				}
			}
			if (add == InvalidHandle())
				std::cout << __FILE__ << ":" << __LINE__ << " bad edges, number of edges " << input.size() << std::endl;
			else
			{
				ret.push_back(add);
				m->SetPrivateMarker(add, mrk);
			}
		}
		ret.RemPrivateMarker(mrk);
		return ret;
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
674
	}
675

Kirill Terekhov's avatar
Kirill Terekhov committed
676
	
677
678
679
	//2d point with comparison operator for line sweep algorithm
	struct Point
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
680
		HandleType node;
Kirill Terekhov's avatar
Kirill Terekhov committed
681
		Storage::real x, y;
Kirill Terekhov's avatar
Kirill Terekhov committed
682
683
684
685
		Point & operator = (Point const & b) { node = b.node; x = b.x; y = b.y; return *this; }
		Point(const Point & b) : node(b.node), x(b.x), y(b.y) {}
		Point(HandleType _node, Storage::real _x, Storage::real _y) : node(_node), x(_x), y(_y) {}
		Point(HandleType _node, Storage::real v[2]) : node(_node), x(v[0]), y(v[1]) {}
686
687
		bool operator <(const Point & b) const
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
688
689
690
			if (y < b.y - ECL_POINT_EPS) return true;
			else if (y > b.y + ECL_POINT_EPS) return false;
			else if (x < b.x - ECL_POINT_EPS) return true;
691
692
693
694
			else return false;
		}
		bool operator ==(const Point & b) const
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
695
			return std::abs(y - b.y) < ECL_POINT_EPS && std::abs(x - b.x) < ECL_POINT_EPS;
696
697
698
		}
		bool operator !=(const Point & b) const
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
699
			return std::abs(y - b.y) > ECL_POINT_EPS || std::abs(x - b.x) > ECL_POINT_EPS;
700
701
		}
	};
702

Kirill Terekhov's avatar
Kirill Terekhov committed
703
	// Used to sort arrays by indices
704
705
706
707
708
709
	class index_comparator
	{
		Storage::integer_array & data;
	public:
		index_comparator(Storage::integer_array & data) : data(data) {}
		index_comparator(const index_comparator & b) : data(b.data) {}
710
711
		index_comparator & operator =(index_comparator const & b) { data = b.data; return *this; }
		bool operator ()(int a, int b) const { return data[a] < data[b]; }
712
	};
Kirill Terekhov's avatar
Kirill Terekhov committed
713
714
715
716
717
718
719
	// Used to sort arrays by depth
	class depth_comparator
	{
		const Storage::real * data;
	public:
		depth_comparator(const Storage::real * data) : data(data) {}
		depth_comparator(const depth_comparator & b) : data(b.data) {}
720
721
		depth_comparator & operator =(depth_comparator const & b) { data = b.data; return *this; }
		bool operator ()(int a, int b) const { return data[a] < data[b]; }
Kirill Terekhov's avatar
Kirill Terekhov committed
722
	};
723

724
725
726
727
	template<typename T>
	int count_duplicates(ElementArray<T> & array)
	{
		Mesh * m = array.GetMeshLink();
728
		MarkerType mrk = m->CreatePrivateMarker();
729
		int dups = 0;
730
		for (int k = 0; k < array.size(); ++k)
731
		{
732
			if (array[k].GetPrivateMarker(mrk))
733
				dups++;
734
			array[k].SetPrivateMarker(mrk);
735
		}
736
737
		array.RemPrivateMarker(mrk);
		m->ReleasePrivateMarker(mrk);
738
739
		return dups;
	}
740

741
742
743
744
	template<typename T>
	void make_unique(ElementArray<T> & array)
	{
		Mesh * m = array.GetMeshLink();
745
		MarkerType mrk = m->CreatePrivateMarker();
746
		int i = 0, j = 0;
747
		while (j < array.size())
748
		{
749
			if (array[j].GetPrivateMarker(mrk))
750
751
752
				++j;
			else
			{
753
				array[j].SetPrivateMarker(mrk);
754
755
756
				array[i++] = array[j++];
			}
		}
757
		array.RemPrivateMarker(mrk);
758
		array.resize(i);
759
		m->ReleasePrivateMarker(mrk);
760
	}
761

Kirill Terekhov's avatar
Kirill Terekhov committed
762
763
	Point make_point(Node n, Tag pnt)
	{
764
		return Point(n->GetHandle(), n->RealArray(pnt).data());
Kirill Terekhov's avatar
Kirill Terekhov committed
765
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
766

767
	std::pair<Storage::real, Storage::real> intersect_pairs(const Point & pabeg, const Point & paend, const Point & pbbeg, const Point & pbend)
Kirill Terekhov's avatar
Kirill Terekhov committed
768
769
	{
		const Storage::real eps = ECL_INTERSECT_EPS;
770
771
		Point pfind(InvalidHandle(), 0, 0);
		Storage::real t1, t2;
Kirill Terekhov's avatar
Kirill Terekhov committed
772
773
		Storage::real div = (pabeg.x - paend.x)*(pbbeg.y - pbend.y) - (pabeg.y - paend.y)*(pbbeg.x - pbend.x);
		if (std::abs(div) < 1.0e-50)
774
			return std::make_pair(-1, -1);
Kirill Terekhov's avatar
Kirill Terekhov committed
775
776
777
		pfind.x = ((pabeg.x*paend.y - pabeg.y*paend.x)*(pbbeg.x - pbend.x) - (pabeg.x - paend.x)*(pbbeg.x*pbend.y - pbbeg.y*pbend.x)) / div;
		pfind.y = ((pabeg.x*paend.y - pabeg.y*paend.x)*(pbbeg.y - pbend.y) - (pabeg.y - paend.y)*(pbbeg.x*pbend.y - pbbeg.y*pbend.x)) / div;
		//optimization! - should belong unit cube
778
		if (pfind.x < 0 || pfind.x > 1 || pfind.y < 0 || pfind.y > 1) return std::make_pair(-1, -1);
Kirill Terekhov's avatar
Kirill Terekhov committed
779
780
781
		if (std::abs(paend.x - pabeg.x) > eps)
		{
			t1 = (pfind.x - pabeg.x) / (paend.x - pabeg.x);
782
			if (t1 < eps || t1 > 1.0 - eps)  { return std::make_pair(-1, -1); }
Kirill Terekhov's avatar
Kirill Terekhov committed
783
784
785
786
		}
		if (std::abs(paend.y - pabeg.y) > eps)
		{
			t1 = (pfind.y - pabeg.y) / (paend.y - pabeg.y);
787
			if (t1 < eps || t1 > 1.0 - eps)  { return std::make_pair(-1, -1); }
Kirill Terekhov's avatar
Kirill Terekhov committed
788
789
790
791
		}
		if (std::abs(pbend.x - pbbeg.x) > eps)
		{
			t2 = (pfind.x - pbbeg.x) / (pbend.x - pbbeg.x);
792
			if (t2 < eps || t2 > 1.0 - eps)  { return std::make_pair(-1, -1); }
Kirill Terekhov's avatar
Kirill Terekhov committed
793
794
795
796
		}
		if (std::abs(pbend.y - pbbeg.y) > eps)
		{
			t2 = (pfind.y - pbbeg.y) / (pbend.y - pbbeg.y);
797
			if (t2 < eps || t2 > 1.0 - eps)  { return std::make_pair(-1, -1); }
Kirill Terekhov's avatar
Kirill Terekhov committed
798
		}
799
		return std::make_pair(t1, t2);
Kirill Terekhov's avatar
Kirill Terekhov committed
800
	}
801

Kirill Terekhov's avatar
Kirill Terekhov committed
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
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


	//returns distance between line if shortest line is within segments, writes position of distance into last argument
	static double SegmentDistance(const double v1[3], const double v2[3], const double v3[3], const double v4[3], double vout[3])
	{
		double v13[3], v21[3], v43[3];
		for (int k = 0; k < 3; ++k)
		{
			v13[k] = v1[k] - v3[k];
			v21[k] = v2[k] - v1[k];
			v43[k] = v4[k] - v3[k];
		}
		double d1321 = 0, d2121 = 0, d4321 = 0, d1343 = 0, d4343 = 0;
		for (int k = 0; k < 3; ++k)
		{
			d1321 += v13[k] * v21[k];
			d2121 += v21[k] * v21[k];
			d4321 += v43[k] * v21[k];
			d1343 += v13[k] * v43[k];
			d4343 += v43[k] * v43[k];
		}
		double mu1 = 0, mu2 = 0;
		/*
		double parallel = 0;
		for (int k = 0; k < 3; ++k)
		parallel += v21[k] * v43[k];
		parallel = fabs(parallel);
		parallel /= sqrt(d2121)*sqrt(d4343);
		if (fabs(parallel - 1.0) > 1.0e-6)
		*/
		if (fabs(d2121*d4343 - d4321*d4321) > 1.0e-6)
		{
			mu1 = (d1343*d4321 - d1321*d4343) / (d2121*d4343 - d4321*d4321);
			mu2 = (d1343 + mu1*d4321) / d4343;
			if (mu1 > 1) mu1 = 1;
			if (mu1 < 0) mu1 = 0;
			if (mu2 > 1) mu2 = 1;
			if (mu2 < 0) mu2 = 0;
		}
		else //parallel lines
		{
			//1d problem
			double la = 0, ra = 0;
			double lb = 0, rb = 0;
			for (int k = 0; k < 3; ++k)
			{
				la += v1[k] * v21[k];
				ra += v2[k] * v21[k];
				lb += v3[k] * v21[k];
				rb += v4[k] * v21[k];
			}
			bool invb = false;
			if (lb > rb)
			{
				std::swap(lb, rb);
				invb = true;
			}

			if (ra+1.0e-7 < lb) // (la,ra) is to the left of (lb,rb), no intersection
			{
				mu1 = 1;
				mu2 = 0;
			}
			else if (la > rb+1.0e-7) // (la,ra) is to the right of (lb,rb), no intersection
			{
				mu1 = 0;
				mu2 = 1;
			}
			else if (lb > la) // (lb,rb) intersects to the right of or contains (la,ra)
			{
				mu2 = 0;
				mu1 = (lb - la) / (ra - la);
			}
			else if (rb < ra) // (lb,rb) intersects to the left of or contains (la,ra)
			{
				mu2 = 1;
				mu1 = (rb - la) / (ra - la);
			}
			else // (lb,rb) contains (la,ra)
			{
				mu1 = 0;
				mu2 = (la - lb) / (rb - lb);
			}
			if (invb) mu2 = 1 - mu2;
		}
		double vs1[3], vs2[3], h = 0;
		for (int k = 0; k < 3; ++k)
		{
			vs1[k] = v1[k] + mu1*(v2[k] - v1[k]);
			vs2[k] = v3[k] + mu2*(v4[k] - v3[k]);
			vout[k] = (vs1[k] + vs2[k])*0.5;
			h += (vs2[k] - vs1[k])*(vs2[k] - vs1[k]);
		}
		return sqrt(h);
	}



	std::pair<bool, Point> intersect_segments(Mesh * m, const Edge & a, const Edge & b, std::set<Point> & intersections, Tag pnt, bool print)
901
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
902
		const Storage::real eps = ECL_INTERSECT_EPS;
903
904
905
906
907
908
909
910
911
912
		Point pfind(InvalidHandle(), 0, 0);
		if (a->getBeg() == b->getBeg() ||
			a->getBeg() == b->getEnd() ||
			a->getEnd() == b->getEnd() ||
			a->getEnd() == b->getBeg())
			return std::make_pair(false, pfind);
		Point pabeg = make_point(a->getBeg(), pnt);
		Point paend = make_point(a->getEnd(), pnt);
		Point pbbeg = make_point(b->getBeg(), pnt);
		Point pbend = make_point(b->getEnd(), pnt);
Kirill Terekhov's avatar
Kirill Terekhov committed
913
		Storage::real div = (pabeg.x - paend.x)*(pbbeg.y - pbend.y) - (pabeg.y - paend.y)*(pbbeg.x - pbend.x);
914
		Storage::real t1, t2;
Kirill Terekhov's avatar
Kirill Terekhov committed
915
		if (std::abs(div) < 1.0e-50)
916
917
		{
			if (print) std::cout << "divisor is zero" << std::endl;
918
			return std::make_pair(false, pfind);
919
920
921
		}
		pfind.x = ((pabeg.x*paend.y - pabeg.y*paend.x)*(pbbeg.x - pbend.x) - (pabeg.x - paend.x)*(pbbeg.x*pbend.y - pbbeg.y*pbend.x)) / div;
		pfind.y = ((pabeg.x*paend.y - pabeg.y*paend.x)*(pbbeg.y - pbend.y) - (pabeg.y - paend.y)*(pbbeg.x*pbend.y - pbbeg.y*pbend.x)) / div;
922
		//optimization uses information that we stay in unit cube
923
924
925
		if (pfind.x < 0 || pfind.x > 1 || pfind.y < 0 || pfind.y > 1)
			return std::make_pair(false, pfind);
		if (print) std::cout << "found (" << pfind.x << ", " << pfind.y << ") for edges " << a->GetHandle() << " and " << b->GetHandle() << std::endl;
926
		//probably some of these tests are redundant
Kirill Terekhov's avatar
Kirill Terekhov committed
927
		if (std::abs(paend.x - pabeg.x) > eps)
928
929
		{
			t1 = (pfind.x - pabeg.x) / (paend.x - pabeg.x);
930
			if (t1 < eps || t1 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false, pfind); }
931
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
932
		if (std::abs(paend.y - pabeg.y) > eps)
933
934
		{
			t1 = (pfind.y - pabeg.y) / (paend.y - pabeg.y);
935
			if (t1 < eps || t1 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false, pfind); }
936
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
937
		if (std::abs(pbend.x - pbbeg.x) > eps)
938
939
		{
			t2 = (pfind.x - pbbeg.x) / (pbend.x - pbbeg.x);
940
			if (t2 < eps || t2 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false, pfind); }
941
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
942
		if (std::abs(pbend.y - pbbeg.y) > eps)
943
944
		{
			t2 = (pfind.y - pbbeg.y) / (pbend.y - pbbeg.y);
945
			if (t2 < eps || t2 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false, pfind); }
946
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
947
		std::set<Point>::iterator search = intersections.find(pfind);
948
		//check whether intersection already exists
949
		if (search != intersections.end()) //there is a node!
Kirill Terekhov's avatar
Kirill Terekhov committed
950
951
		{
			pfind = *search;
952
953
954
955
956
			if (pfind.node == a->getBeg()->GetHandle() ||
				pfind.node == b->getBeg()->GetHandle() ||
				pfind.node == a->getEnd()->GetHandle() ||
				pfind.node == b->getEnd()->GetHandle())
				return std::make_pair(false, pfind);
Kirill Terekhov's avatar
Kirill Terekhov committed
957
		}
958
959
		else //no node, create one
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
960
			Storage::real find[3];
961
			//restore coordinate
Kirill Terekhov's avatar
Kirill Terekhov committed
962
963
			SegmentDistance(a->getBeg().Coords().data(), a->getEnd().Coords().data(), b->getBeg().Coords().data(), b->getEnd().Coords().data(), find);

Kirill Terekhov's avatar
Kirill Terekhov committed
964
965
			pfind.node = m->CreateNode(find)->GetHandle();
			if (print) std::cout << "intersection accepted (" << find[0] << "," << find[1] << "," << find[2] << ") t1 " << t1 << " t2 " << t2 << " new node " << pfind.node << std::endl;
966
			Storage::real_array _pfind = m->RealArray(pfind.node, pnt);
967
968
			_pfind[0] = pfind.x;
			_pfind[1] = pfind.y;
Kirill Terekhov's avatar
Kirill Terekhov committed
969
			intersections.insert(pfind);
970
		}
971
		return std::make_pair(true, pfind);
972
	}
973

Kirill Terekhov's avatar
Kirill Terekhov committed
974
975
976
977
	void split_edge(Mesh * m, Node I, Edge a, ElementArray<Edge> & splitted_a, std::vector<Tag> & transfer, bool print)
	{
		std::vector< std::vector<char> > copy(transfer.size());
		//memorize data
978
		for (int k = 0; k < transfer.size(); ++k) //tags
Kirill Terekhov's avatar
Kirill Terekhov committed
979
980
		{
			int size = a->GetDataSize(transfer[k]);
981
982
			copy[k].resize(transfer[k].GetBytesSize()*size);
			if (!copy.empty()) a->GetData(transfer[k], 0, size, &copy[k][0]);
Kirill Terekhov's avatar
Kirill Terekhov committed
983
		}
984
985
986
987

		if (print) std::cout << "split a " << a->GetHandle() << " " << a->getBeg()->GetHandle() << " <-> " << a->getEnd()->GetHandle() << ":-> ";
		splitted_a = Edge::SplitEdge(a, ElementArray<Node>(m, 1, I->GetHandle()), 0);
		if (print) std::cout << splitted_a[0]->GetHandle() << " " << splitted_a[0]->getBeg()->GetHandle() << " <-> " << splitted_a[0]->getEnd()->GetHandle() << " and " << splitted_a[1]->GetHandle() << " " << splitted_a[1]->getBeg()->GetHandle() << " <-> " << splitted_a[1]->getEnd()->GetHandle() << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
988
		//duplicate data
989
		for (int k = 0; k < transfer.size(); ++k)
Kirill Terekhov's avatar
Kirill Terekhov committed
990
991
		{
			int size = (int)copy[k].size() / transfer[k].GetBytesSize();
992
			if (size) for (int l = 0; l < 2; ++l) //two parts
Kirill Terekhov's avatar
Kirill Terekhov committed
993
			{
994
995
				splitted_a[l].SetDataSize(transfer[k], size);
				splitted_a[l].SetData(transfer[k], 0, size, &copy[k][0]);
Kirill Terekhov's avatar
Kirill Terekhov committed
996
997
998
			}
		}
	}
999

1000
	void split_edges(Mesh * m, Node I, Edge a, Edge b, ElementArray<Edge> & splitted_a, ElementArray<Edge> & splitted_b, std::vector<Tag> & transfer, bool print)
1001
	{
1002
1003
		split_edge(m, I, a, splitted_a, transfer, print);
		split_edge(m, I, b, splitted_b, transfer, print);
1004
	}
1005

Kirill Terekhov's avatar
Kirill Terekhov committed
1006
1007
	
	void intersect_naive(Mesh * m, ElementArray<Edge> & segments, ElementArray<Node> & nodes, std::vector<Tag> & transfer, Tag pnt, bool print)
1008
	{
1009
		//Tag pnt = m->CreateTag("PROJ_PNT"+m->GetLocalProcessorRank(),DATA_REAL,NODE,NODE,2);
Kirill Terekhov's avatar
Kirill Terekhov committed
1010
		std::set<Point> intersections;
1011
		std::vector<HandleType> initials(segments.size() * 2);
1012
		MarkerType initial = m->CreatePrivateMarker();
1013
1014
		for (int k = 0; k < (int)segments.size(); ++k)
		{
1015
1016
			initials[k * 2 + 0] = segments[k]->getBeg()->GetHandle();
			initials[k * 2 + 1] = segments[k]->getEnd()->GetHandle();
1017
1018
			segments[k]->getBeg()->SetPrivateMarker(initial);
			segments[k]->getEnd()->SetPrivateMarker(initial);
1019
1020
			intersections.insert(make_point(segments[k]->getBeg(), pnt));
			intersections.insert(make_point(segments[k]->getEnd(), pnt));
1021
		}
1022
		for (int i = 0; i < (int)segments.size(); ++i)
1023
		{
1024
			for (int j = i + 1; j < (int)segments.size(); ++j)
1025
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1026
				std::pair<bool, Point> I = intersect_segments(m, segments[i], segments[j], intersections, pnt, print);
1027
				if (I.first)
1028
1029
				{
					ElementArray<Edge> splitted_a, splitted_b;
Kirill Terekhov's avatar
Kirill Terekhov committed
1030
1031
1032
					split_edges(m, Node(m, I.second.node), segments[i], segments[j], splitted_a, splitted_b, transfer, print);
					segments[i] = splitted_a[0];
					segments[j] = splitted_b[0];
1033
1034
1035
1036
1037
					segments.push_back(splitted_a[1]);
					segments.push_back(splitted_b[1]);
				}
			}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
1038
		//nodes.clear();
1039
		for (std::set<Point>::iterator it = intersections.begin(); it != intersections.end(); ++it)
1040
		{
1041
			if (!m->GetPrivateMarker(it->node, initial))
Kirill Terekhov's avatar
Kirill Terekhov committed
1042
				nodes.push_back(it->node);
1043
		}
1044
1045
		for (int k = 0; k < (int)initials.size(); ++k)
			m->RemPrivateMarker(initials[k], initial);
1046
		m->ReleasePrivateMarker(initial);
1047
1048
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
1049
	void block_number_union(Element n, ElementArray<Element> & adj, Tag block, Tag write)
1050
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
1051
		Storage::integer_array bn = n->IntegerArray(write);
1052
		if (!adj.empty())
1053
		{
1054
1055
1056
			dynarray<int, 64> uni, tmp;
			//std::vector<int> uni, tmp;
			for (ElementArray<Edge>::size_type k = 0; k < adj.size(); ++k)
1057
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1058
				Storage::integer_array be = adj[k]->IntegerArray(block);
1059
1060
				tmp.resize(uni.size() + be.size());
				tmp.resize(std::set_union(uni.begin(), uni.end(), be.begin(), be.end(), tmp.begin()) - tmp.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1061
				uni.swap(tmp);
1062
			}
1063
			bn.replace(bn.begin(), bn.end(), uni.begin(), uni.end());
1064
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
1065
		else bn.clear();
1066
	}
1067

Kirill Terekhov's avatar
Kirill Terekhov committed
1068
1069
1070
	void block_number_union_merge(Element n, ElementArray<Element> & adj, Tag block, Tag write)
	{
		Storage::integer_array bn = n->IntegerArray(write);
1071
		if (!adj.empty())
1072
		{
1073
1074
1075
			dynarray<int, 64> uni(bn.begin(), bn.end()), tmp(uni.size());
			//std::vector<int> uni(bn.begin(), bn.end()), tmp(uni.size());
			for (ElementArray<Edge>::size_type k = 0; k < adj.size(); ++k)
1076
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1077
				Storage::integer_array be = adj[k]->IntegerArray(block);
1078
1079
				tmp.resize(uni.size() + be.size());
				tmp.resize(std::set_union(uni.begin(), uni.end(), be.begin(), be.end(), tmp.begin()) - tmp.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1080
				uni.swap(tmp);
1081
			}
1082
			bn.replace(bn.begin(), bn.end(), uni.begin(), uni.end());
1083
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
1084
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
1085

1086
	void block_number_intersection(Element n, const ElementArray<Element> & adj, Tag block, Tag write)
Kirill Terekhov's avatar
Kirill Terekhov committed
1087
1088
	{
		Storage::integer_array bn = n->IntegerArray(write);
1089
		if (!adj.empty())
Kirill Terekhov's avatar
Kirill Terekhov committed
1090
1091
		{
			Storage::integer_array be = adj[0]->IntegerArray(block);
1092
1093
1094
			dynarray<int, 64> inter(be.begin(), be.end()), tmp(inter.size());
			//std::vector<int> inter(be.begin(), be.end()), tmp(inter.size());
			for (ElementArray<Edge>::size_type k = 1; k < adj.size(); ++k)
Kirill Terekhov's avatar
Kirill Terekhov committed
1095
1096
			{
				be = adj[k]->IntegerArray(block);
1097
				tmp.resize(std::set_intersection(inter.begin(), inter.end(), be.begin(), be.end(), tmp.begin()) - tmp.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1098
1099
				inter.swap(tmp);
			}
1100
			bn.replace(bn.begin(), bn.end(), inter.begin(), inter.end());
Kirill Terekhov's avatar
Kirill Terekhov committed
1101
1102
1103
1104
1105
1106
		}
		else bn.clear();
		//Storage::integer_array bn = n->IntegerArray(block);
		//bn.replace(bn.begin(),bn.end(),inter.begin(),inter.end());
	}

1107
1108

	struct mult_rec
Kirill Terekhov's avatar
Kirill Terekhov committed
1109
	{
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
		std::string name;
		double mult;
		int bijk[6];
	};

	struct editnnc_rec
	{
		int bijk[6];
		double mult;
	};

	struct fault_rec
	{
		std::string name;
		int dir;
		int bijk[6];
	};

	struct multflt_rec
	{
		std::string name;
		double mult;
	};
1133

Kirill Terekhov's avatar
Kirill Terekhov committed
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
	struct prod_rate
	{
		double orat;
		double wrat;
		double grat;
		double lrat;
		double resv;
	};

	struct inje_rate
	{
		double rate;
		double resv;
	};

	struct wconprodinje_rec
	{
		std::string type;
		std::string ctrl;
		bool open;
		union
		{
			prod_rate prat;
			inje_rate irat;
			double rats[5];
		} urats;
		double bhp;
	};
	
	struct compdat_rec
	{
		int bijk[4];
		bool open;
		int satnum;
		double WI, skin, perm, dfac,rw,r0;
		char dir;
	};
	struct welspecs_rec
	{
		int bij[2];
		double depth;
		std::string group;
		std::string phase;
	};
	typedef std::vector<compdat_rec> compdat_recs;
	struct well_rec
	{
		welspecs_rec wspec;
		std::map<int,wconprodinje_rec> wpi;
		std::map<int,compdat_recs> compdat;
	};
	//typedef std::map<std::string, compdat_recs > compdat_wells;
	//typedef std::map<std::string, welspecs_rec > welspecs_wells;
	typedef std::map<std::string, well_rec > wells_data;
	struct schedule
	{
		char tagnum;
		double value;
		HandleType element;
	};
	typedef std::vector<schedule> schedules;



1198
1199
	void Mesh::LoadECL(std::string File)
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
1200
		char have_perm = 0;
1201
		std::cout << std::scientific;
Kirill Terekhov's avatar
Kirill Terekhov committed
1202
		bool perform_splitting = false;
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
1203
1204
		bool project_perm = false;
		int split_degenerate = 0;
1205
1206
1207
		bool check_topology = false;
		int verbosity = 0;
		for (INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); ++k)
Kirill Terekhov's avatar
Kirill Terekhov committed
1208
		{
1209
			if (file_options[k].first == "ECL_SPLIT_GLUED")
Kirill Terekhov's avatar
Kirill Terekhov committed
1210
			{
1211
				if (file_options[k].second == "TRUE")
Kirill Terekhov's avatar
Kirill Terekhov committed
1212
1213
1214
1215
					perform_splitting = true;
				else
					perform_splitting = false;
			}
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
			if (file_options[k].first == "ECL_DEGENERATE") //TODO4
			{
				if (file_options[k].second == "TRUE")
					split_degenerate = 1;
				else if (file_options[k].second == "TRANM")
					split_degenerate = 2;
				else
					split_degenerate = 0;
			}
			if (file_options[k].first == "ECL_TOPOLOGY")
Kirill Terekhov's avatar
Kirill Terekhov committed
1226
			{
1227
1228
1229
1230
1231
1232
1233
1234
				if (file_options[k].second == "TRUE")
					check_topology = true;
				else
					check_topology = false;
			}
			if (file_options[k].first == "ECL_PROJECT_PERM")
			{
				if (file_options[k].second == "TRUE")
Kirill Terekhov's avatar
Kirill Terekhov committed
1235
1236
1237
1238
					project_perm = true;
				else
					project_perm = false;
			}
1239
1240
1241
1242
1243
1244
1245
1246
1247
			if (file_options[k].first == "VERBOSITY")
			{
				verbosity = atoi(file_options[k].second.c_str());
				if (verbosity < 0 || verbosity > 2)
				{
					printf("%s:%d Unknown verbosity option: %s\n", __FILE__, __LINE__, file_options[k].second.c_str());
					verbosity = 1;
				}
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
1248
1249
1250
1251
		}
		std::vector<HandleType> old_nodes(NumberOfNodes());
		{
			unsigned qq = 0;
1252
			for (Mesh::iteratorNode it = BeginNode(); it != EndNode(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
1253
1254
				old_nodes[qq++] = *it;
		}
1255
1256
		if (!old_nodes.empty())
			std::sort(old_nodes.begin(), old_nodes.end(), CentroidComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1257

1258
1259
		FILE * f = fopen(File.c_str(), "r");
		if (f == NULL)
Kirill Terekhov's avatar
Kirill Terekhov committed
1260
1261
1262
1263
		{
			std::cout << __FILE__ << ":" << __LINE__ << " cannot open file " << File << std::endl;
			throw BadFileName;
		}
1264
1265
1266
1267
1268
1269
1270
		double tt;
		if (verbosity)
		{
			tt = Timer();
			std::cout << "Started loading " << File << std::endl;
		}
		std::vector< std::pair< std::pair<FILE *, std::string>, int> > fs(1, std::make_pair(std::make_pair(f, File), 0));
Kirill Terekhov's avatar
Kirill Terekhov committed
1271
		char readline[2048], readlines[2048], *p, *pend, rec[2048];
1272
		int text_end, text_start, state = ECL_NONE, state_from = ECL_NONE, nchars;
Kirill Terekhov's avatar
Kirill Terekhov committed
1273
1274
1275
1276
1277
1278
1279
		int waitlines = 0;
		int have_dimens = 0, totread, downread, numrecs, offset;
		int gtype = ECL_GTYPE_NONE;
		int argtype = ECL_VAR_NONE;
		int radial = ECL_GTYPE_NONE;
		Storage::real * read_arrayf = NULL;
		Storage::integer * read_arrayi = NULL;
1280
		Storage::integer dims[3], mapaxis[6] = { 0, 1, 0, 0, 1, 0 };
Kirill Terekhov's avatar
Kirill Terekhov committed
1281
		Storage::real inrad = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1282
		std::vector<Storage::real> xyz, perm, poro, tops, zcorn, ntg, pressure, pbub, swat, soil, sgas, thconr;
1283
1284
1285
1286
1287
1288
1289