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

#include "inmost.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
11
#include "../Mesh/incident_matrix.hpp"
Kirill Terekhov's avatar
Kirill Terekhov committed
12
#include <string>
Kirill Terekhov's avatar
Kirill Terekhov committed
13
#include <ctime>
Kirill Terekhov's avatar
Kirill Terekhov committed
14
15
#if defined(USE_MESH)

Kirill Terekhov's avatar
Kirill Terekhov committed
16
17
18
19
20
21
22
23
24
25
// 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
26
//todo
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
27
28
// 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
29
// 2. fix fast intersect() algorithm, check against intersect_naive()
30
31
32
33
// 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
34
// 5. (ok) populate keywords data, poro, perm, actnum, satnum, etc...
35
36
37
// 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
38
// 6. local grid refinement
39
40
// 6.1 radial local grid refinement
// 6.2 nested local grid refinement
41
// 7. (ok) omp-parallel loading
42
// 8. (ok) mpi-parallel loading
43
// 9. (ok) do not convert arrays xyz and zcorn
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
44
//10. read wells into sets
Kirill Terekhov's avatar
Kirill Terekhov committed
45
//10.1 (ok) read compdat
Kirill Terekhov's avatar
Kirill Terekhov committed
46
47
48
49
//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
50
51
//11. optimize:
//11.1 detect that pillars have no faults, skip intersection in this case
52
53
54
55
//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
56
57
//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
58
59
60
61
//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
62

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

Kirill Terekhov's avatar
Kirill Terekhov committed
65
66
67
68
69
70

//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
71
72
73
//this controls what is recorded into tag CTRL of well sets
#define ECL_WCTRL_RATE 0
#define ECL_WCTRL_BHP 1
74
#define ECL_WCTRL_RESV 2
Kirill Terekhov's avatar
Kirill Terekhov committed
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

//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
94
#define ECL_WTAG_BHP 7
Kirill Terekhov's avatar
Kirill Terekhov committed
95

96
//eclipse states
Kirill Terekhov's avatar
Kirill Terekhov committed
97
98
99
100
101
102
103
104
105
106
107
108
#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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
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
161
162
163
164
165
166
167
#define ECL_PERMXY 11
#define ECL_PERMXZ 12
#define ECL_PERMYZ 13
#define ECL_PORO 14
#define ECL_MAPAXIS 15
#define ECL_INRAD 16
#define ECL_COORDS 17
#define ECL_ZCORN 18
#define ECL_ACTNUM 19
#define ECL_SATNUM 20
#define ECL_THCONR 21
#define ECL_MULTIPLY 22
#define ECL_MULTIPLY_MUL 23
#define ECL_MULTIPLY_BLK 24
#define ECL_EDITNNC 25
#define ECL_EDITNNC_BLK 26
#define ECL_EDITNNC_MUL 27
#define ECL_NTG 28
#define ECL_FAULTS 29
#define ECL_FAULTS_BLK 30
#define ECL_FAULTS_DIR 31
#define ECL_MULTFLT 32
#define ECL_MULTFLT_MUL 33
#define ECL_PRESSURE 34
#define ECL_SWAT 35
#define ECL_SOIL 36
#define ECL_SGAS 37
#define ECL_PBUB 38
#define ECL_EQLNUM 39
#define ECL_ROCKNUM 40
#define ECL_PVTNUM 41
#define ECL_COMPDAT 42
#define ECL_COMPDAT_BLK 43
#define ECL_COMPDAT_OPEN 44
#define ECL_COMPDAT_SATNUM 45
#define ECL_COMPDAT_TRANS 46
#define ECL_COMPDAT_BORE 47
#define ECL_COMPDAT_PERM 48
#define ECL_COMPDAT_SKIN 49
#define ECL_COMPDAT_DFAC 50
#define ECL_COMPDAT_DIR 51
#define ECL_COMPDAT_RAD 52
#define ECL_DATES 53
#define ECL_DATES_MON 54
#define ECL_DATES_YEAR 55
#define ECL_DATES_HRS 56
#define ECL_WELSPECS 57
#define ECL_WELSPECS_GROUP 58
#define ECL_WELSPECS_I 59
#define ECL_WELSPECS_J 60
#define ECL_WELSPECS_Z 61
#define ECL_WELSPECS_PHASE 62
#define ECL_TSTEP 63
#define ECL_WCONPRODINJE 64
#define ECL_WCONPRODINJE_TYPE 65
#define ECL_WCONPRODINJE_OPEN 66
#define ECL_WCONPRODINJE_CTRL 67
#define ECL_WCONPRODINJE_RATES 68
#define ECL_WCONPRODINJE_BHP 69
168

Kirill Terekhov's avatar
Kirill Terekhov committed
169
170
171
172
173
174
175
176
177
178
179
180

#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

181
#define ECL_IJK_DATA(i,j,k) (i + ((j)+(k)*dims[1])*dims[0])
182
183
#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)
184
185
186
187
188

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

Kirill Terekhov's avatar
Kirill Terekhov committed
189
190
191
#define HAVE_PERM_X 1
#define HAVE_PERM_Y 2
#define HAVE_PERM_Z 4
192
193
194
#define HAVE_PERM_XY 8
#define HAVE_PERM_XZ 16
#define HAVE_PERM_YZ 32
Kirill Terekhov's avatar
Kirill Terekhov committed
195

196

197

Kirill Terekhov's avatar
Kirill Terekhov committed
198
199
namespace INMOST
{
200

Kirill Terekhov's avatar
Kirill Terekhov committed
201

202
	void GetDXYZ(Element e, INMOST_DATA_REAL_TYPE & dx, INMOST_DATA_REAL_TYPE & dy, INMOST_DATA_REAL_TYPE & dz)
Kirill Terekhov's avatar
Kirill Terekhov committed
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
270
271
272
273
274
275
276
277
278
279
	{
		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;
	}
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
	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;
302
			for (size_t k = 0; k < n_faces.size(); ++k) if (n_faces[k].GetMarker(frac_markers)) num_frac_faces++;
303
304
305
			if (num_frac_faces)
			{
				ElementArray<Cell> n_cells = it->getCells(); //get all cells of the node
306
				for (size_t k = 0; k < n_cells.size(); ++k) //assign local enumeration to the cells
307
					n_cells[k].IntegerDF(indexes) = static_cast<Storage::integer>(k);
308
				for (size_t k = 0; k < n_faces.size(); ++k) //stich together node's numbers corresponding to cells if no fracture separates them
309
310
311
312
313
314
315
				{
					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)
						{
316
							for (size_t q = 0; q < n_cells.size(); q++) if (n_cells[q].IntegerDF(indexes) == fi)
317
318
319
320
321
								n_cells[q].IntegerDF(indexes) = bi;
						}
					}
				}
				dynarray<int, 64> nums(n_cells.size()); //store all numbers
322
				for (size_t k = 0; k < n_cells.size(); ++k) nums[k] = n_cells[k].IntegerDF(indexes);
323
324
325
326
				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
				{
327
					for (size_t k = 0; k < nums.size(); k++)
328
329
330
331
					{
						Node image;
						Storage::real xyz[3] = { 0, 0, 0 }, cntc[3] = { 0, 0, 0 };
						int n_node_cells = 0;
332
						for (size_t q = 0; q < n_cells.size(); q++)
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
						{
							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);

351
						for (size_t q = 0; q < n_cells.size(); q++)
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
						{
							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
432
				for (size_t q = 0; q < nodes.size(); ++q) if (nodes[q].GetMarker(frac_markers))
433
434
435
436
437
438
439
440
441
				{
					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
442
				for (size_t k = 0; k < edges.size(); k++)
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
				{
					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
464
				for (size_t q = 0; q < nodes.size(); ++q) if (nodes[q].GetMarker(frac_markers))
465
466
467
468
469
470
471
472
473
				{
					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
474
				for (size_t k = 0; k < edges.size(); k++)
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
521
522
523
524
525
526
527
528
529
530
				{
					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
531
	//replace n* witn n stars
532
	void UnStar(const char * input, char * output)
Kirill Terekhov's avatar
Kirill Terekhov committed
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
561
562
563
564
565
566
567
568
569
570
	{
		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
		}
	}

571
	std::string STR8(const char * input)
572
573
574
575
576
577
578
579
	{
		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;
	}

580
	std::string UnQoute(std::string input)
581
582
583
584
585
586
587
	{
		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;
	}
588
	std::string ToUpper(std::string input)
589
590
591
592
593
594
	{
		for (size_t k = 0; k < input.size(); ++k)
			input[k] = toupper(input[k]);
		return input;
	}

595
	int GetDir(std::string dir)
596
597
598
	{
		if (dir == "X-" || dir == "I-")
			return 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
599
		else if (dir == "X" || dir == "I" || dir == "X+" || dir == "I+" || dir == "X " || dir == "I ")
600
601
602
			return 1;
		else if (dir == "Y-" || dir == "J-")
			return 2;
Kirill Terekhov's avatar
Kirill Terekhov committed
603
		else if (dir == "Y" || dir == "J" || dir == "Y+" || dir == "J+" || dir == "Y " || dir == "J ")
604
605
606
			return 3;
		else if (dir == "Z-" || dir == "K-")
			return 4;
Kirill Terekhov's avatar
Kirill Terekhov committed
607
		else if (dir == "Z" || dir == "K" || dir == "Z+" || dir == "K+" || dir == "Z " || dir == "K ")
608
609
610
			return 5;
		return -1;
	}
611
	std::string GetFolder(std::string file)
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
612
613
	{
		size_t found = file.find_last_of("/\\");
614
		if (found == std::string::npos)
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
615
			return "";
616
617
618
		else return file.substr(0, found);
	}

619
	ElementArray<Edge> OrderEdges(const ElementArray<Edge> & input, MarkerType mrk)
620
621
622
623
	{
		Mesh * m = const_cast<Mesh *>(input.GetMeshLink());
		ElementArray<Edge> ret(m);
		ret.reserve(input.size());
624
		HandleType last = InvalidHandle();
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
674
675
676
677
678
679
680
681
682
683
		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
684
	}
685

Kirill Terekhov's avatar
Kirill Terekhov committed
686
	
687
688
689
	//2d point with comparison operator for line sweep algorithm
	struct Point
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
690
		HandleType node;
Kirill Terekhov's avatar
Kirill Terekhov committed
691
		Storage::real x, y;
Kirill Terekhov's avatar
Kirill Terekhov committed
692
693
694
695
		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]) {}
696
697
		bool operator <(const Point & b) const
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
698
699
700
			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;
701
702
703
704
			else return false;
		}
		bool operator ==(const Point & b) const
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
705
			return std::abs(y - b.y) < ECL_POINT_EPS && std::abs(x - b.x) < ECL_POINT_EPS;
706
707
708
		}
		bool operator !=(const Point & b) const
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
709
			return std::abs(y - b.y) > ECL_POINT_EPS || std::abs(x - b.x) > ECL_POINT_EPS;
710
711
		}
	};
712

Kirill Terekhov's avatar
Kirill Terekhov committed
713
	// Used to sort arrays by indices
714
715
716
717
718
719
	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) {}
720
721
		index_comparator & operator =(index_comparator const & b) { data = b.data; return *this; }
		bool operator ()(int a, int b) const { return data[a] < data[b]; }
722
	};
Kirill Terekhov's avatar
Kirill Terekhov committed
723
724
725
726
727
728
729
	// 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) {}
730
731
		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
732
	};
733

734
735
736
737
	template<typename T>
	int count_duplicates(ElementArray<T> & array)
	{
		Mesh * m = array.GetMeshLink();
738
		MarkerType mrk = m->CreatePrivateMarker();
739
		int dups = 0;
740
		for (typename ElementArray<T>::size_type k = 0; k < array.size(); ++k)
741
		{
742
			if (array[k].GetPrivateMarker(mrk))
743
				dups++;
744
			array[k].SetPrivateMarker(mrk);
745
		}
746
747
		array.RemPrivateMarker(mrk);
		m->ReleasePrivateMarker(mrk);
748
749
		return dups;
	}
750

751
752
753
754
	template<typename T>
	void make_unique(ElementArray<T> & array)
	{
		Mesh * m = array.GetMeshLink();
755
		MarkerType mrk = m->CreatePrivateMarker();
756
		typename ElementArray<T>::size_type i = 0, j = 0;
757
		while (j < array.size())
758
		{
759
			if (array[j].GetPrivateMarker(mrk))
760
761
762
				++j;
			else
			{
763
				array[j].SetPrivateMarker(mrk);
764
765
766
				array[i++] = array[j++];
			}
		}
767
		array.RemPrivateMarker(mrk);
768
		array.resize(i);
769
		m->ReleasePrivateMarker(mrk);
770
	}
771

Kirill Terekhov's avatar
Kirill Terekhov committed
772
773
	Point make_point(Node n, Tag pnt)
	{
774
		return Point(n->GetHandle(), n->RealArray(pnt).data());
Kirill Terekhov's avatar
Kirill Terekhov committed
775
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
776

777
	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
778
779
	{
		const Storage::real eps = ECL_INTERSECT_EPS;
780
		Point pfind(InvalidHandle(), 0, 0);
781
		Storage::real t1 = -1, t2 = -1;
Kirill Terekhov's avatar
Kirill Terekhov committed
782
783
		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)
784
			return std::make_pair(-1, -1);
Kirill Terekhov's avatar
Kirill Terekhov committed
785
786
787
		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
788
		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
789
790
791
		if (std::abs(paend.x - pabeg.x) > eps)
		{
			t1 = (pfind.x - pabeg.x) / (paend.x - pabeg.x);
792
			if (t1 < eps || t1 > 1.0 - eps)  { return std::make_pair(-1, -1); }
Kirill Terekhov's avatar
Kirill Terekhov committed
793
794
795
796
		}
		if (std::abs(paend.y - pabeg.y) > eps)
		{
			t1 = (pfind.y - pabeg.y) / (paend.y - pabeg.y);
797
			if (t1 < eps || t1 > 1.0 - eps)  { return std::make_pair(-1, -1); }
Kirill Terekhov's avatar
Kirill Terekhov committed
798
799
800
801
		}
		if (std::abs(pbend.x - pbbeg.x) > eps)
		{
			t2 = (pfind.x - pbbeg.x) / (pbend.x - pbbeg.x);
802
			if (t2 < eps || t2 > 1.0 - eps)  { return std::make_pair(-1, -1); }
Kirill Terekhov's avatar
Kirill Terekhov committed
803
804
805
806
		}
		if (std::abs(pbend.y - pbbeg.y) > eps)
		{
			t2 = (pfind.y - pbbeg.y) / (pbend.y - pbbeg.y);
807
			if (t2 < eps || t2 > 1.0 - eps)  { return std::make_pair(-1, -1); }
Kirill Terekhov's avatar
Kirill Terekhov committed
808
		}
809
		return std::make_pair(t1, t2);
Kirill Terekhov's avatar
Kirill Terekhov committed
810
	}
811

Kirill Terekhov's avatar
Kirill Terekhov committed
812
813
814


	//returns distance between line if shortest line is within segments, writes position of distance into last argument
815
	static INMOST_DATA_REAL_TYPE SegmentDistance(const INMOST_DATA_REAL_TYPE v1[3], const INMOST_DATA_REAL_TYPE v2[3], const INMOST_DATA_REAL_TYPE v3[3], const INMOST_DATA_REAL_TYPE v4[3], INMOST_DATA_REAL_TYPE vout[3])
Kirill Terekhov's avatar
Kirill Terekhov committed
816
	{
817
		INMOST_DATA_REAL_TYPE v13[3], v21[3], v43[3];
Kirill Terekhov's avatar
Kirill Terekhov committed
818
819
820
821
822
823
		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];
		}
824
		INMOST_DATA_REAL_TYPE d1321 = 0, d2121 = 0, d4321 = 0, d1343 = 0, d4343 = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
825
826
827
828
829
830
831
832
		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];
		}
833
		INMOST_DATA_REAL_TYPE mu1 = 0, mu2 = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
		/*
		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
854
855
			INMOST_DATA_REAL_TYPE la = 0, ra = 0;
			INMOST_DATA_REAL_TYPE lb = 0, rb = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
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
			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;
		}
897
		INMOST_DATA_REAL_TYPE vs1[3], vs2[3], h = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
898
899
900
901
902
903
904
905
906
907
908
909
910
		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)
911
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
912
		const Storage::real eps = ECL_INTERSECT_EPS;
913
914
915
916
917
918
919
920
921
922
		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
923
		Storage::real div = (pabeg.x - paend.x)*(pbbeg.y - pbend.y) - (pabeg.y - paend.y)*(pbbeg.x - pbend.x);
924
		Storage::real t1, t2;
Kirill Terekhov's avatar
Kirill Terekhov committed
925
		if (std::abs(div) < 1.0e-50)
926
927
		{
			if (print) std::cout << "divisor is zero" << std::endl;
928
			return std::make_pair(false, pfind);
929
930
931
		}
		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;
932
		//optimization uses information that we stay in unit cube
933
934
935
		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;
936
		//probably some of these tests are redundant
Kirill Terekhov's avatar
Kirill Terekhov committed
937
		if (std::abs(paend.x - pabeg.x) > eps)
938
939
		{
			t1 = (pfind.x - pabeg.x) / (paend.x - pabeg.x);
940
			if (t1 < eps || t1 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false, pfind); }
941
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
942
		if (std::abs(paend.y - pabeg.y) > eps)
943
944
		{
			t1 = (pfind.y - pabeg.y) / (paend.y - pabeg.y);
945
			if (t1 < eps || t1 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false, pfind); }
946
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
947
		if (std::abs(pbend.x - pbbeg.x) > eps)
948
949
		{
			t2 = (pfind.x - pbbeg.x) / (pbend.x - pbbeg.x);
950
			if (t2 < eps || t2 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false, pfind); }
951
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
952
		if (std::abs(pbend.y - pbbeg.y) > eps)
953
954
		{
			t2 = (pfind.y - pbbeg.y) / (pbend.y - pbbeg.y);
955
			if (t2 < eps || t2 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false, pfind); }
956
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
957
		std::set<Point>::iterator search = intersections.find(pfind);
958
		//check whether intersection already exists
959
		if (search != intersections.end()) //there is a node!
Kirill Terekhov's avatar
Kirill Terekhov committed
960
961
		{
			pfind = *search;
962
963
964
965
966
			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
967
		}
968
969
		else //no node, create one
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
970
			Storage::real find[3];
971
			//restore coordinate
Kirill Terekhov's avatar
Kirill Terekhov committed
972
973
			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
974
975
			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;
976
			Storage::real_array _pfind = m->RealArray(pfind.node, pnt);
977
978
			_pfind[0] = pfind.x;
			_pfind[1] = pfind.y;
Kirill Terekhov's avatar
Kirill Terekhov committed
979
			intersections.insert(pfind);
980
		}
981
		return std::make_pair(true, pfind);
982
	}
983

Kirill Terekhov's avatar
Kirill Terekhov committed
984
985
986
987
	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
988
		for (size_t k = 0; k < transfer.size(); ++k) //tags
Kirill Terekhov's avatar
Kirill Terekhov committed
989
		{
990
			INMOST_DATA_ENUM_TYPE size = a->GetDataSize(transfer[k]);
991
			copy[k].resize((size_t)transfer[k].GetBytesSize()*size);
992
			if (!copy.empty()) a->GetData(transfer[k], 0, size, &copy[k][0]);
Kirill Terekhov's avatar
Kirill Terekhov committed
993
		}
994
995
996
997

		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
998
		//duplicate data
999
		for (size_t k = 0; k < transfer.size(); ++k)
Kirill Terekhov's avatar
Kirill Terekhov committed
1000
		{
1001
			INMOST_DATA_ENUM_TYPE size = static_cast<INMOST_DATA_ENUM_TYPE>(copy[k].size() / transfer[k].GetBytesSize());
1002
			if (size) for (int l = 0; l < 2; ++l) //two parts
Kirill Terekhov's avatar
Kirill Terekhov committed
1003
			{
1004
1005
				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
1006
1007
1008
			}
		}
	}
1009

1010
	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)
1011
	{
1012
1013
		split_edge(m, I, a, splitted_a, transfer, print);
		split_edge(m, I, b, splitted_b, transfer, print);
1014
	}
1015

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

Kirill Terekhov's avatar
Kirill Terekhov committed
1059
	void block_number_union(Element n, ElementArray<Element> & adj, Tag block, Tag write)
1060
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
1061
		Storage::integer_array bn = n->IntegerArray(write);
1062
		if (!adj.empty())
1063
		{
1064
1065
1066
			dynarray<int, 64> uni, tmp;
			//std::vector<int> uni, tmp;
			for (ElementArray<Edge>::size_type k = 0; k < adj.size(); ++k)
1067
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1068
				Storage::integer_array be = adj[k]->IntegerArray(block);
1069
1070
				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
1071
				uni.swap(tmp);
1072
			}
1073
			bn.replace(bn.begin(), bn.end(), uni.begin(), uni.end());
1074
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
1075
		else bn.clear();
1076
	}
1077

Kirill Terekhov's avatar
Kirill Terekhov committed
1078
1079
1080
	void block_number_union_merge(Element n, ElementArray<Element> & adj, Tag block, Tag write)
	{
		Storage::integer_array bn = n->IntegerArray(write);
1081
		if (!adj.empty())
1082
		{
1083
1084
1085
			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)
1086
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1087
				Storage::integer_array be = adj[k]->IntegerArray(block);
1088
1089
				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
1090
				uni.swap(tmp);
1091
			}
1092
			bn.replace(bn.begin(), bn.end(), uni.begin(), uni.end());
1093
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
1094
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
1095

1096
	void block_number_intersection(Element n, const ElementArray<Element> & adj, Tag block, Tag write)
Kirill Terekhov's avatar
Kirill Terekhov committed
1097
1098
	{
		Storage::integer_array bn = n->IntegerArray(write);
1099
		if (!adj.empty())
Kirill Terekhov's avatar
Kirill Terekhov committed
1100
1101
		{
			Storage::integer_array be = adj[0]->IntegerArray(block);
1102
1103
1104
			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
1105
1106
			{
				be = adj[k]->IntegerArray(block);
1107
				tmp.resize(std::set_intersection(inter.begin(), inter.end(), be.begin(), be.end(), tmp.begin()) - tmp.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1108
1109
				inter.swap(tmp);
			}
1110
			bn.replace(bn.begin(), bn.end(), inter.begin(), inter.end());
Kirill Terekhov's avatar
Kirill Terekhov committed
1111
1112
1113
1114
1115
1116
		}
		else bn.clear();
		//Storage::integer_array bn = n->IntegerArray(block);
		//bn.replace(bn.begin(),bn.end(),inter.begin(),inter.end());
	}

1117
1118

	struct mult_rec
Kirill Terekhov's avatar
Kirill Terekhov committed
1119
	{
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
		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;
	};
1143

Kirill Terekhov's avatar
Kirill Terekhov committed
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
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
	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;



1208
1209
	void Mesh::LoadECL(std::string File)
	{
Kirill Terekhov's avatar
sync    
Kirill Terekhov committed
1210
		bool parallel_read = true;
Kirill Terekhov's avatar
Kirill Terekhov committed
1211
		char have_perm = 0;
1212
		std::cout << std::scientific;
1213
		int perform_splitting = 0;
Kirill Terekhov's avatar
Updates