mesh_ecl_file.cpp 191 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
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>
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
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
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
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
Kirill Terekhov committed
44
//10. read wells into sets
Kirill Terekhov's avatar
Kirill Terekhov committed
45
//10.1 (ok) read compdat
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
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

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
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
95

96
//eclipse states
Kirill Terekhov's avatar
Kirill Terekhov committed
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
#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
114
#define ECL_ACTNUM 16
Kirill Terekhov's avatar
Kirill Terekhov committed
115
#define ECL_SATNUM 17
116
#define ECL_THCONR 18
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
#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
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
#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
165

Kirill Terekhov's avatar
Kirill Terekhov committed
166 167 168 169 170 171 172 173 174 175 176 177

#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

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

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

Kirill Terekhov's avatar
Kirill Terekhov committed
186 187 188
#define HAVE_PERM_X 1
#define HAVE_PERM_Y 2
#define HAVE_PERM_Z 4
Kirill Terekhov's avatar
Kirill Terekhov committed
189

190

191

Kirill Terekhov's avatar
Kirill Terekhov committed
192 193
namespace INMOST
{
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 270 271 272 273

	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;
	}
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 521 522 523 524
	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");
	}

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 561 562 563 564
	//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
		}
	}

565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592
	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
593
		else if (dir == "X" || dir == "I" || dir == "X+" || dir == "I+" || dir == "X " || dir == "I ")
594 595 596
			return 1;
		else if (dir == "Y-" || dir == "J-")
			return 2;
Kirill Terekhov's avatar
Kirill Terekhov committed
597
		else if (dir == "Y" || dir == "J" || dir == "Y+" || dir == "J+" || dir == "Y " || dir == "J ")
598 599 600
			return 3;
		else if (dir == "Z-" || dir == "K-")
			return 4;
Kirill Terekhov's avatar
Kirill Terekhov committed
601
		else if (dir == "Z" || dir == "K" || dir == "Z+" || dir == "K+" || dir == "Z " || dir == "K ")
602 603 604
			return 5;
		return -1;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
605 606 607
	static std::string GetFolder(std::string file)
	{
		size_t found = file.find_last_of("/\\");
608
		if (found == std::string::npos)
Kirill Terekhov's avatar
Kirill Terekhov committed
609
			return "";
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 674 675 676 677
		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
Kirill Terekhov committed
678
	}
679

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

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

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

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

Kirill Terekhov's avatar
Kirill Terekhov committed
766 767
	Point make_point(Node n, Tag pnt)
	{
768
		return Point(n->GetHandle(), n->RealArray(pnt).data());
Kirill Terekhov's avatar
Kirill Terekhov committed
769
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
770

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


	//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)
905
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
906
		const Storage::real eps = ECL_INTERSECT_EPS;
907 908 909 910 911 912 913 914 915 916
		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
917
		Storage::real div = (pabeg.x - paend.x)*(pbbeg.y - pbend.y) - (pabeg.y - paend.y)*(pbbeg.x - pbend.x);
918
		Storage::real t1, t2;
Kirill Terekhov's avatar
Kirill Terekhov committed
919
		if (std::abs(div) < 1.0e-50)
920 921
		{
			if (print) std::cout << "divisor is zero" << std::endl;
922
			return std::make_pair(false, pfind);
923 924 925
		}
		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;
926
		//optimization uses information that we stay in unit cube
927 928 929
		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;
930
		//probably some of these tests are redundant
Kirill Terekhov's avatar
Kirill Terekhov committed
931
		if (std::abs(paend.x - pabeg.x) > eps)
932 933
		{
			t1 = (pfind.x - pabeg.x) / (paend.x - pabeg.x);
934
			if (t1 < eps || t1 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false, pfind); }
935
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
936
		if (std::abs(paend.y - pabeg.y) > eps)
937 938
		{
			t1 = (pfind.y - pabeg.y) / (paend.y - pabeg.y);
939
			if (t1 < eps || t1 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false, pfind); }
940
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
941
		if (std::abs(pbend.x - pbbeg.x) > eps)
942 943
		{
			t2 = (pfind.x - pbbeg.x) / (pbend.x - pbbeg.x);
944
			if (t2 < eps || t2 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false, pfind); }
945
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
946
		if (std::abs(pbend.y - pbbeg.y) > eps)
947 948
		{
			t2 = (pfind.y - pbbeg.y) / (pbend.y - pbbeg.y);
949
			if (t2 < eps || t2 > 1.0 - eps)  { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false, pfind); }
950
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
951
		std::set<Point>::iterator search = intersections.find(pfind);
952
		//check whether intersection already exists
953
		if (search != intersections.end()) //there is a node!
Kirill Terekhov's avatar
Kirill Terekhov committed
954 955
		{
			pfind = *search;
956 957 958 959 960
			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
961
		}
962 963
		else //no node, create one
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
964
			Storage::real find[3];
965
			//restore coordinate
966 967
			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
968 969
			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;
970
			Storage::real_array _pfind = m->RealArray(pfind.node, pnt);
971 972
			_pfind[0] = pfind.x;
			_pfind[1] = pfind.y;
Kirill Terekhov's avatar
Kirill Terekhov committed
973
			intersections.insert(pfind);
974
		}
975
		return std::make_pair(true, pfind);
976
	}
977

Kirill Terekhov's avatar
Kirill Terekhov committed
978 979 980 981
	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
982
		for (int k = 0; k < transfer.size(); ++k) //tags
Kirill Terekhov's avatar
Kirill Terekhov committed
983 984
		{
			int size = a->GetDataSize(transfer[k]);
985 986
			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
987
		}
988 989 990 991

		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
992
		//duplicate data
993
		for (int k = 0; k < transfer.size(); ++k)
Kirill Terekhov's avatar
Kirill Terekhov committed
994 995
		{
			int size = (int)copy[k].size() / transfer[k].GetBytesSize();
996
			if (size) for (int l = 0; l < 2; ++l) //two parts
Kirill Terekhov's avatar
Kirill Terekhov committed
997
			{
998 999
				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
1000 1001 1002
			}
		}
	}
1003

1004
	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)
1005
	{
1006 1007
		split_edge(m, I, a, splitted_a, transfer, print);
		split_edge(m, I, b, splitted_b, transfer, print);
1008
	}
1009

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

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

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

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

1111 1112

	struct mult_rec
Kirill Terekhov's avatar
Kirill Terekhov committed
1113
	{
1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136
		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;
	};
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 1198 1199 1200 1201
	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;



1202 1203
	void Mesh::LoadECL(std::string File)
	{
Kirill Terekhov's avatar
sync  
Kirill Terekhov committed
1204
		bool parallel_read = true;
Kirill Terekhov's avatar
Kirill Terekhov committed
1205
		char have_perm = 0;
1206
		std::cout << std::scientific;
1207
		int perform_splitting = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1208 1209
		bool project_perm = false;
		int split_degenerate = 0;
1210 1211 1212
		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
1213
		{
1214
			if (file_options[k].first == "ECL_SPLIT_GLUED")
Kirill Terekhov's avatar
Kirill Terekhov committed
1215
			{
1216 1217 1218
				if (file_options[k].second == "ALL")
					perform_splitting = 2;
				else if (file_options[k].second == "TRUE")
1219
					perform_splitting = 1;
Kirill Terekhov's avatar
Kirill Terekhov committed
1220
				else
1221
					perform_splitting = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1222
			}
1223 1224 1225 1226 1227 1228 1229 1230 1231 1232
			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
1233
			{
1234 1235 1236 1237 1238 1239 1240 1241
				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
1242 1243 1244 1245
					project_perm = true;
				else
					project_perm = false;
			}
1246 1247 1248 1249 1250 1251 1252 1253 1254
			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
1255 1256 1257 1258
		}
		std::vector<HandleType> old_nodes(NumberOfNodes());
		{
			unsigned qq = 0;
1259
			for (Mesh::iteratorNode it = BeginNode(); it != EndNode(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
1260 1261
				old_nodes[qq++] = *it;
		}
1262 1263
		if (!old_nodes.empty())
			std::sort(old_nodes.begin(), old_nodes.end(), CentroidComparator(this));
Kirill Terekhov's avatar
Kirill Terekhov committed
1264

1265 1266
		FILE * f = fopen(File.c_str(), "r");
		if (f == NULL)
Kirill Terekhov's avatar
Kirill Terekhov committed
1267 1268 1269 1270
		{
			std::cout << __FILE__ << ":" << __LINE__ << " cannot open file " << File << std::endl;
			throw BadFileName;
		}
1271 1272 1273 1274 1275 1276 1277
		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));
1278
		char readline[2048], readlines[2048], *p, *pend, rec[2048], pupper[2048];
1279
		int text_end, text_start, state = ECL_NONE, state_from = ECL_NONE, state_incl = ECL_NONE, nchars;
Kirill Terekhov's avatar
Kirill Terekhov committed
1280 1281 1282 1283 1284 1285 1286
		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;
1287
		Storage::integer dims[3], mapaxis[6] = { 0, 1, 0, 0, 1, 0 };
Kirill Terekhov's avatar
Kirill Terekhov committed
1288
		Storage::real inrad = 0;
1289
		std::vector<Storage::real> xyz, perm, poro, tops, zcorn, ntg, pressure, pbub, swat, soil, sgas, thconr;
1290 1291 1292 1293 1294 1295 1296 1297 1298
		std::vector<Storage::integer> actnum, satnum, eqlnum, pvtnum, rocknum;
		std::vector<mult_rec> multiply;
		mult_rec mult_cur;
		std::vector<editnnc_rec> editnnc;
		editnnc_rec editnnc_cur;
		std::vector<fault_rec> faults;
		fault_rec fault_cur;
		std::vector<multflt_rec> multflt;
		multflt_rec multflt_cur;
1299 1300 1301 1302 1303 1304 1305 1306
		time_t start, last;
		bool have_start = false;
		struct tm date_cur;
		bool read_start;
		std::vector<double> tsteps(1,0.0); //skip first tstep
		std::pair<std::string,compdat_rec> compdat_cur;
		std::pair<std::string,welspecs_rec> welspecs_cur;
		std::pair<std::string,wconprodinje_rec> wconprodinje_cur;
1307 1308
		bool multiply_trans = false;
		bool multiply_props = false;
1309 1310 1311 1312 1313
		bool wconprod = false; //WCONPROD or WCONINJE
		bool wellspec = false;// WELLSPEC keyword instead of WELLSPECS
		//compdat_wells compdat; 
		//welspecs_wells welspecs;
		wells_data wells_sched;
1314
		while (!fs.empty())
Kirill Terekhov's avatar
Kirill Terekhov committed
1315
		{
1316
			while (fgets(readline, 2048, fs.back().first.first) != NULL)
Kirill Terekhov's avatar
Kirill Terekhov committed
1317 1318 1319
			{
				fs.back().second++; //line number
				{
1320
					if (readline[strlen(readline) - 1] == '\n') readline[strlen(readline) - 1] = '\0';
1321 1322
					UnStar(readline,readlines); //remove n*
					strcpy(readline,readlines); //copy back
Kirill Terekhov's avatar
Kirill Terekhov committed
1323
					text_end = static_cast<int>(strlen(readline));
1324 1325 1326 1327
					for (text_start = 0; isspace(readline[text_start]) && text_start < text_end; text_start++);
					if (text_start == text_end) continue;
					for (text_end = text_end - 1; isspace(readline[text_end]) && text_end > text_start; text_end--);
					readline[text_end + 1] = '\0';
Kirill Terekhov's avatar
Kirill Terekhov committed
1328 1329
					p = readline + text_start;
					pend = readline + text_end + 1;
1330 1331
					strcpy(pupper,p); //upper case for keyword comparison
					for (char * q = pupper; q < pupper + (pend-p); q++) *q = toupper(*q);
Kirill Terekhov's avatar
Kirill Terekhov committed
1332
				}
1333 1334 1335
				if (p[0] == '-' && p[1] == '-') continue; //skip comment
				if (waitlines) { waitlines--; continue; } //skip meaningful lines
				switch (state)
Kirill Terekhov's avatar
Kirill Terekhov committed
1336 1337
				{
				case ECL_NONE:
1338
					if (!ECLSTRCMP(pupper, "END")) //end of data - don't read beyond
Kirill Terekhov's avatar
Kirill Terekhov committed
1339 1340 1341
					{
						goto ecl_exit_loop;
					}
1342 1343
					else if (!ECLSTRCMP(pupper, "INCLUDE")) state = ECL_INCLUDE;
					else if (!ECLSTRCMP(pupper, "BOX"))
Kirill Terekhov's avatar
Kirill Terekhov committed
1344
					{