mesh_ecl_file.cpp 191 KB
 Kirill Terekhov committed Nov 23, 2015 1 ``````#ifdef _MSC_VER //kill some warnings `````` igor committed Mar 28, 2018 2 ``````#ifndef _CRT_SECURE_NO_WARNINGS `````` Kirill Terekhov committed Nov 23, 2015 3 ``````#define _CRT_SECURE_NO_WARNINGS `````` igor committed Mar 28, 2018 4 5 ``````#endif #ifndef _SCL_SECURE_NO_WARNINGS `````` Kirill Terekhov committed Jun 21, 2017 6 ``````#define _SCL_SECURE_NO_WARNINGS `````` Kirill Terekhov committed Nov 23, 2015 7 ``````#endif `````` igor committed Mar 28, 2018 8 ``````#endif `````` Kirill Terekhov committed Nov 23, 2015 9 10 `````` #include "inmost.h" `````` Kirill Terekhov committed Jul 09, 2016 11 ``````#include "../Mesh/incident_matrix.hpp" `````` Kirill Terekhov committed Oct 04, 2016 12 ``````#include `````` Kirill Terekhov committed Jun 21, 2017 13 ``````#include `````` Kirill Terekhov committed Nov 23, 2015 14 15 ``````#if defined(USE_MESH) `````` Kirill Terekhov committed Jul 09, 2016 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 committed Jul 18, 2016 26 ``````//todo `````` Kirill Terekhov committed Jul 27, 2016 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 committed Jul 18, 2016 29 ``````// 2. fix fast intersect() algorithm, check against intersect_naive() `````` Kirill Terekhov committed Jun 15, 2017 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 committed Jul 27, 2016 34 ``````// 5. (ok) populate keywords data, poro, perm, actnum, satnum, etc... `````` Kirill Terekhov committed Jun 15, 2017 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 committed Jul 27, 2016 38 ``````// 6. local grid refinement `````` Kirill Terekhov committed Jul 29, 2016 39 40 ``````// 6.1 radial local grid refinement // 6.2 nested local grid refinement `````` Kirill Terekhov committed Jul 28, 2016 41 ``````// 7. (ok) omp-parallel loading `````` Kirill Terekhov committed Jun 15, 2017 42 ``````// 8. (ok) mpi-parallel loading `````` Kirill Terekhov committed Jul 29, 2016 43 ``````// 9. (ok) do not convert arrays xyz and zcorn `````` Kirill Terekhov committed Jul 27, 2016 44 ``````//10. read wells into sets `````` Kirill Terekhov committed Nov 03, 2016 45 ``````//10.1 (ok) read compdat `````` Kirill Terekhov committed Jun 21, 2017 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 `````` Kirill Terekhov committed Jul 29, 2016 50 51 ``````//11. optimize: //11.1 detect that pillars have no faults, skip intersection in this case `````` Kirill Terekhov committed Jun 15, 2017 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 committed Jun 21, 2017 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 `````` Kirill Terekhov committed Jun 15, 2017 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 committed Jul 18, 2016 62 `````` `````` Kirill Terekhov committed Jun 04, 2017 63 ``````#define ECLSTRCMP(x,y) strncmp(STR8(x).c_str(),STR8(y).c_str(),8) `````` Kirill Terekhov committed Nov 23, 2015 64 `````` `````` Kirill Terekhov committed Aug 12, 2016 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 committed Jun 21, 2017 71 72 73 ``````//this controls what is recorded into tag CTRL of well sets #define ECL_WCTRL_RATE 0 #define ECL_WCTRL_BHP 1 `````` Kirill Terekhov committed Jun 23, 2017 74 ``````#define ECL_WCTRL_RESV 2 `````` Kirill Terekhov committed Jun 21, 2017 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 `````` Kirill Terekhov committed Jun 23, 2017 94 ``````#define ECL_WTAG_BHP 7 `````` Kirill Terekhov committed Jun 21, 2017 95 `````` `````` Kirill Terekhov committed Jun 15, 2017 96 ``````//eclipse states `````` Kirill Terekhov committed Nov 23, 2015 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 `````` Kirill Terekhov committed Jul 18, 2016 114 ``````#define ECL_ACTNUM 16 `````` Kirill Terekhov committed Jul 27, 2016 115 ``````#define ECL_SATNUM 17 `````` Kirill Terekhov committed Jun 21, 2017 116 ``````#define ECL_THCONR 18 `````` Kirill Terekhov committed Jun 15, 2017 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 `````` Kirill Terekhov committed Jun 21, 2017 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 `````` Kirill Terekhov committed Jun 15, 2017 165 `````` `````` Kirill Terekhov committed Nov 23, 2015 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 `````` Kirill Terekhov committed Jul 08, 2016 178 ``````#define ECL_IJK_DATA(i,j,k) (i + ((j)+(k)*dims[1])*dims[0]) `````` Kirill Terekhov committed Jul 29, 2016 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) `````` Kirill Terekhov committed Jul 08, 2016 181 182 183 184 185 `````` //line sweep events #define SEG_START 0 #define SEG_END 1 `````` Kirill Terekhov committed May 28, 2017 186 187 188 ``````#define HAVE_PERM_X 1 #define HAVE_PERM_Y 2 #define HAVE_PERM_Z 4 `````` Kirill Terekhov committed Nov 23, 2015 189 `````` `````` Kirill Terekhov committed Jun 15, 2017 190 `````` `````` Kirill Terekhov committed Jun 04, 2017 191 `````` `````` Kirill Terekhov committed Nov 23, 2015 192 193 ``````namespace INMOST { `````` Kirill Terekhov committed Jun 15, 2017 194 `````` `````` Kirill Terekhov committed Jun 21, 2017 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 nodes = e->getNodes(); for (ElementArray::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; } `````` Kirill Terekhov committed Jun 15, 2017 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 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 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 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(n_node_cells); xyz[1] /= static_cast(n_node_cells); xyz[2] /= static_cast(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 edges = it->getEdges(); //go through edges of current face ElementArray newedges(m); Storage::reference_array images = it->BackCell()->ReferenceArray(cell2node); //mark me and all my edges that they may contain images for (ElementArray::iterator jt = edges.begin(); jt != edges.end(); ++jt) { if (jt->getBeg()->GetMarker(frac_markers) || jt->getEnd()->GetMarker(frac_markers)) { ElementArray enodes(m, 2); if (jt->getBeg()->GetMarker(frac_markers)) //search image of node between my elements { for (int k = 0; k < static_cast(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(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 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 fracfaces(m); ElementArray nodes = it->getNodes(), nodesb(m, nodes.size()), nodesf(m, nodes.size()); ElementArray 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 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 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 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 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 change(m); ElementArray faces = it->getFaces(); for (ElementArray::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 committed Jun 21, 2017 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 } } `````` Kirill Terekhov committed Jun 15, 2017 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 committed Aug 05, 2017 593 `````` else if (dir == "X" || dir == "I" || dir == "X+" || dir == "I+" || dir == "X " || dir == "I ") `````` Kirill Terekhov committed Jun 15, 2017 594 595 596 `````` return 1; else if (dir == "Y-" || dir == "J-") return 2; `````` Kirill Terekhov committed Aug 05, 2017 597 `````` else if (dir == "Y" || dir == "J" || dir == "Y+" || dir == "J+" || dir == "Y " || dir == "J ") `````` Kirill Terekhov committed Jun 15, 2017 598 599 600 `````` return 3; else if (dir == "Z-" || dir == "K-") return 4; `````` Kirill Terekhov committed Aug 05, 2017 601 `````` else if (dir == "Z" || dir == "K" || dir == "Z+" || dir == "K+" || dir == "Z " || dir == "K ") `````` Kirill Terekhov committed Jun 15, 2017 602 603 604 `````` return 5; return -1; } `````` Kirill Terekhov committed Jul 27, 2016 605 606 607 `````` static std::string GetFolder(std::string file) { size_t found = file.find_last_of("/\\"); `````` Kirill Terekhov committed Jun 15, 2017 608 `````` if (found == std::string::npos) `````` Kirill Terekhov committed Jul 27, 2016 609 `````` return ""; `````` Kirill Terekhov committed Jun 15, 2017 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 OrderEdges(const ElementArray & input, MarkerType mrk) { Mesh * m = const_cast(input.GetMeshLink()); ElementArray 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 committed Jul 27, 2016 678 `````` } `````` Kirill Terekhov committed Jun 15, 2017 679 `````` `````` Kirill Terekhov committed Jun 21, 2017 680 `````` `````` Kirill Terekhov committed Jul 08, 2016 681 682 683 `````` //2d point with comparison operator for line sweep algorithm struct Point { `````` Kirill Terekhov committed Aug 12, 2016 684 `````` HandleType node; `````` Kirill Terekhov committed Nov 04, 2016 685 `````` Storage::real x, y; `````` Kirill Terekhov committed Aug 12, 2016 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]) {} `````` Kirill Terekhov committed Jul 08, 2016 690 691 `````` bool operator <(const Point & b) const { `````` Kirill Terekhov committed Aug 12, 2016 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; `````` Kirill Terekhov committed Jul 08, 2016 695 696 697 698 `````` else return false; } bool operator ==(const Point & b) const { `````` Kirill Terekhov committed Aug 12, 2016 699 `````` return std::abs(y - b.y) < ECL_POINT_EPS && std::abs(x - b.x) < ECL_POINT_EPS; `````` Kirill Terekhov committed Jul 08, 2016 700 701 702 `````` } bool operator !=(const Point & b) const { `````` Kirill Terekhov committed Aug 12, 2016 703 `````` return std::abs(y - b.y) > ECL_POINT_EPS || std::abs(x - b.x) > ECL_POINT_EPS; `````` Kirill Terekhov committed Jul 29, 2016 704 705 `````` } }; `````` Kirill Terekhov committed Jun 15, 2017 706 `````` `````` Kirill Terekhov committed Aug 12, 2016 707 `````` // Used to sort arrays by indices `````` Kirill Terekhov committed Jul 18, 2016 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) {} `````` Kirill Terekhov committed Jun 15, 2017 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]; } `````` Kirill Terekhov committed Jul 18, 2016 716 `````` }; `````` Kirill Terekhov committed Aug 12, 2016 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) {} `````` Kirill Terekhov committed Jun 15, 2017 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 committed Aug 12, 2016 726 `````` }; `````` Kirill Terekhov committed Jun 15, 2017 727 `````` `````` Kirill Terekhov committed Jul 18, 2016 728 729 730 731 `````` template int count_duplicates(ElementArray & array) { Mesh * m = array.GetMeshLink(); `````` Kirill Terekhov committed Jul 28, 2016 732 `````` MarkerType mrk = m->CreatePrivateMarker(); `````` Kirill Terekhov committed Jul 18, 2016 733 `````` int dups = 0; `````` Kirill Terekhov committed Jun 15, 2017 734 `````` for (int k = 0; k < array.size(); ++k) `````` Kirill Terekhov committed Jul 18, 2016 735 `````` { `````` Kirill Terekhov committed Jun 15, 2017 736 `````` if (array[k].GetPrivateMarker(mrk)) `````` Kirill Terekhov committed Jul 18, 2016 737 `````` dups++; `````` Kirill Terekhov committed Jul 28, 2016 738 `````` array[k].SetPrivateMarker(mrk); `````` Kirill Terekhov committed Jul 18, 2016 739 `````` } `````` Kirill Terekhov committed Jul 28, 2016 740 741 `````` array.RemPrivateMarker(mrk); m->ReleasePrivateMarker(mrk); `````` Kirill Terekhov committed Jul 18, 2016 742 743 `````` return dups; } `````` Kirill Terekhov committed Jun 15, 2017 744 `````` `````` Kirill Terekhov committed Jul 18, 2016 745 746 747 748 `````` template void make_unique(ElementArray & array) { Mesh * m = array.GetMeshLink(); `````` Kirill Terekhov committed Jul 28, 2016 749 `````` MarkerType mrk = m->CreatePrivateMarker(); `````` Kirill Terekhov committed Jul 18, 2016 750 `````` int i = 0, j = 0; `````` Kirill Terekhov committed Jun 15, 2017 751 `````` while (j < array.size()) `````` Kirill Terekhov committed Jul 18, 2016 752 `````` { `````` Kirill Terekhov committed Jun 15, 2017 753 `````` if (array[j].GetPrivateMarker(mrk)) `````` Kirill Terekhov committed Jul 18, 2016 754 755 756 `````` ++j; else { `````` Kirill Terekhov committed Jul 28, 2016 757 `````` array[j].SetPrivateMarker(mrk); `````` Kirill Terekhov committed Jul 18, 2016 758 759 760 `````` array[i++] = array[j++]; } } `````` Kirill Terekhov committed Jul 28, 2016 761 `````` array.RemPrivateMarker(mrk); `````` Kirill Terekhov committed Jul 18, 2016 762 `````` array.resize(i); `````` Kirill Terekhov committed Jul 28, 2016 763 `````` m->ReleasePrivateMarker(mrk); `````` Kirill Terekhov committed Jul 18, 2016 764 `````` } `````` Kirill Terekhov committed Jun 15, 2017 765 `````` `````` Kirill Terekhov committed Aug 12, 2016 766 767 `````` Point make_point(Node n, Tag pnt) { `````` Kirill Terekhov committed Jun 15, 2017 768 `````` return Point(n->GetHandle(), n->RealArray(pnt).data()); `````` Kirill Terekhov committed Aug 12, 2016 769 `````` } `````` Kirill Terekhov committed Nov 04, 2016 770 `````` `````` Kirill Terekhov committed Jun 15, 2017 771 `````` std::pair intersect_pairs(const Point & pabeg, const Point & paend, const Point & pbbeg, const Point & pbend) `````` Kirill Terekhov committed Nov 04, 2016 772 773 `````` { const Storage::real eps = ECL_INTERSECT_EPS; `````` Kirill Terekhov committed Jun 15, 2017 774 775 `````` Point pfind(InvalidHandle(), 0, 0); Storage::real t1, t2; `````` Kirill Terekhov committed Nov 04, 2016 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) `````` Kirill Terekhov committed Jun 15, 2017 778 `````` return std::make_pair(-1, -1); `````` Kirill Terekhov committed Nov 04, 2016 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 `````` Kirill Terekhov committed Jun 15, 2017 782 `````` if (pfind.x < 0 || pfind.x > 1 || pfind.y < 0 || pfind.y > 1) return std::make_pair(-1, -1); `````` Kirill Terekhov committed Nov 04, 2016 783 784 785 `````` if (std::abs(paend.x - pabeg.x) > eps) { t1 = (pfind.x - pabeg.x) / (paend.x - pabeg.x); `````` Kirill Terekhov committed Jun 15, 2017 786 `````` if (t1 < eps || t1 > 1.0 - eps) { return std::make_pair(-1, -1); } `````` Kirill Terekhov committed Nov 04, 2016 787 788 789 790 `````` } if (std::abs(paend.y - pabeg.y) > eps) { t1 = (pfind.y - pabeg.y) / (paend.y - pabeg.y); `````` Kirill Terekhov committed Jun 15, 2017 791 `````` if (t1 < eps || t1 > 1.0 - eps) { return std::make_pair(-1, -1); } `````` Kirill Terekhov committed Nov 04, 2016 792 793 794 795 `````` } if (std::abs(pbend.x - pbbeg.x) > eps) { t2 = (pfind.x - pbbeg.x) / (pbend.x - pbbeg.x); `````` Kirill Terekhov committed Jun 15, 2017 796 `````` if (t2 < eps || t2 > 1.0 - eps) { return std::make_pair(-1, -1); } `````` Kirill Terekhov committed Nov 04, 2016 797 798 799 800 `````` } if (std::abs(pbend.y - pbbeg.y) > eps) { t2 = (pfind.y - pbbeg.y) / (pbend.y - pbbeg.y); `````` Kirill Terekhov committed Jun 15, 2017 801 `````` if (t2 < eps || t2 > 1.0 - eps) { return std::make_pair(-1, -1); } `````` Kirill Terekhov committed Nov 04, 2016 802 `````` } `````` Kirill Terekhov committed Jun 15, 2017 803 `````` return std::make_pair(t1, t2); `````` Kirill Terekhov committed Nov 04, 2016 804 `````` } `````` Kirill Terekhov committed Jun 15, 2017 805 `````` `````` Kirill Terekhov committed Jun 21, 2017 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 intersect_segments(Mesh * m, const Edge & a, const Edge & b, std::set & intersections, Tag pnt, bool print) `````` Kirill Terekhov committed Jul 08, 2016 905 `````` { `````` Kirill Terekhov committed Aug 12, 2016 906 `````` const Storage::real eps = ECL_INTERSECT_EPS; `````` Kirill Terekhov committed Jun 15, 2017 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 committed Nov 04, 2016 917 `````` Storage::real div = (pabeg.x - paend.x)*(pbbeg.y - pbend.y) - (pabeg.y - paend.y)*(pbbeg.x - pbend.x); `````` Kirill Terekhov committed Jun 15, 2017 918 `````` Storage::real t1, t2; `````` Kirill Terekhov committed Aug 12, 2016 919 `````` if (std::abs(div) < 1.0e-50) `````` Kirill Terekhov committed Jul 18, 2016 920 921 `````` { if (print) std::cout << "divisor is zero" << std::endl; `````` Kirill Terekhov committed Jun 15, 2017 922 `````` return std::make_pair(false, pfind); `````` Kirill Terekhov committed Jul 18, 2016 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; `````` Kirill Terekhov committed Jul 28, 2016 926 `````` //optimization uses information that we stay in unit cube `````` Kirill Terekhov committed Jun 15, 2017 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; `````` Kirill Terekhov committed Jul 18, 2016 930 `````` //probably some of these tests are redundant `````` Kirill Terekhov committed Aug 12, 2016 931 `````` if (std::abs(paend.x - pabeg.x) > eps) `````` Kirill Terekhov committed Jul 18, 2016 932 933 `````` { t1 = (pfind.x - pabeg.x) / (paend.x - pabeg.x); `````` Kirill Terekhov committed Jun 15, 2017 934 `````` if (t1 < eps || t1 > 1.0 - eps) { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false, pfind); } `````` Kirill Terekhov committed Jul 18, 2016 935 `````` } `````` Kirill Terekhov committed Aug 12, 2016 936 `````` if (std::abs(paend.y - pabeg.y) > eps) `````` Kirill Terekhov committed Jul 18, 2016 937 938 `````` { t1 = (pfind.y - pabeg.y) / (paend.y - pabeg.y); `````` Kirill Terekhov committed Jun 15, 2017 939 `````` if (t1 < eps || t1 > 1.0 - eps) { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false, pfind); } `````` Kirill Terekhov committed Jul 18, 2016 940 `````` } `````` Kirill Terekhov committed Aug 12, 2016 941 `````` if (std::abs(pbend.x - pbbeg.x) > eps) `````` Kirill Terekhov committed Jul 18, 2016 942 943 `````` { t2 = (pfind.x - pbbeg.x) / (pbend.x - pbbeg.x); `````` Kirill Terekhov committed Jun 15, 2017 944 `````` if (t2 < eps || t2 > 1.0 - eps) { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false, pfind); } `````` Kirill Terekhov committed Jul 18, 2016 945 `````` } `````` Kirill Terekhov committed Aug 12, 2016 946 `````` if (std::abs(pbend.y - pbbeg.y) > eps) `````` Kirill Terekhov committed Jul 18, 2016 947 948 `````` { t2 = (pfind.y - pbbeg.y) / (pbend.y - pbbeg.y); `````` Kirill Terekhov committed Jun 15, 2017 949 `````` if (t2 < eps || t2 > 1.0 - eps) { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false, pfind); } `````` Kirill Terekhov committed Jul 18, 2016 950 `````` } `````` Kirill Terekhov committed Aug 12, 2016 951 `````` std::set::iterator search = intersections.find(pfind); `````` Kirill Terekhov committed Jul 18, 2016 952 `````` //check whether intersection already exists `````` Kirill Terekhov committed Jun 15, 2017 953 `````` if (search != intersections.end()) //there is a node! `````` Kirill Terekhov committed Aug 12, 2016 954 955 `````` { pfind = *search; `````` Kirill Terekhov committed Jun 15, 2017 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 committed Aug 12, 2016 961 `````` } `````` Kirill Terekhov committed Jul 18, 2016 962 963 `````` else //no node, create one { `````` Kirill Terekhov committed Aug 12, 2016 964 `````` Storage::real find[3]; `````` Kirill Terekhov committed Jul 29, 2016 965 `````` //restore coordinate `````` Kirill Terekhov committed Jun 21, 2017 966 967 `````` SegmentDistance(a->getBeg().Coords().data(), a->getEnd().Coords().data(), b->getBeg().Coords().data(), b->getEnd().Coords().data(), find); `````` Kirill Terekhov committed Aug 12, 2016 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; `````` Kirill Terekhov committed Jun 15, 2017 970 `````` Storage::real_array _pfind = m->RealArray(pfind.node, pnt); `````` Kirill Terekhov committed Jul 28, 2016 971 972 `````` _pfind[0] = pfind.x; _pfind[1] = pfind.y; `````` Kirill Terekhov committed Aug 12, 2016 973 `````` intersections.insert(pfind); `````` Kirill Terekhov committed Jul 18, 2016 974 `````` } `````` Kirill Terekhov committed Jun 15, 2017 975 `````` return std::make_pair(true, pfind); `````` Kirill Terekhov committed Jul 18, 2016 976 `````` } `````` Kirill Terekhov committed Jul 08, 2016 977 `````` `````` Kirill Terekhov committed Nov 03, 2016 978 979 980 981 `````` void split_edge(Mesh * m, Node I, Edge a, ElementArray & splitted_a, std::vector & transfer, bool print) { std::vector< std::vector > copy(transfer.size()); //memorize data `````` Kirill Terekhov committed Jun 15, 2017 982 `````` for (int k = 0; k < transfer.size(); ++k) //tags `````` Kirill Terekhov committed Nov 03, 2016 983 984 `````` { int size = a->GetDataSize(transfer[k]); `````` Kirill Terekhov committed Jun 15, 2017 985 986 `````` copy[k].resize(transfer[k].GetBytesSize()*size); if (!copy.empty()) a->GetData(transfer[k], 0, size, ©[k][0]); `````` Kirill Terekhov committed Nov 03, 2016 987 `````` } `````` Kirill Terekhov committed Jun 15, 2017 988 989 990 991 `````` if (print) std::cout << "split a " << a->GetHandle() << " " << a->getBeg()->GetHandle() << " <-> " << a->getEnd()->GetHandle() << ":-> "; splitted_a = Edge::SplitEdge(a, ElementArray(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 committed Nov 03, 2016 992 `````` //duplicate data `````` Kirill Terekhov committed Jun 15, 2017 993 `````` for (int k = 0; k < transfer.size(); ++k) `````` Kirill Terekhov committed Nov 03, 2016 994 995 `````` { int size = (int)copy[k].size() / transfer[k].GetBytesSize(); `````` Kirill Terekhov committed Jun 15, 2017 996 `````` if (size) for (int l = 0; l < 2; ++l) //two parts `````` Kirill Terekhov committed Nov 03, 2016 997 `````` { `````` Kirill Terekhov committed Jun 15, 2017 998 999 `````` splitted_a[l].SetDataSize(transfer[k], size); splitted_a[l].SetData(transfer[k], 0, size, ©[k][0]); `````` Kirill Terekhov committed Nov 03, 2016 1000 1001 1002 `````` } } } `````` Kirill Terekhov committed Jul 18, 2016 1003 `````` `````` Kirill Terekhov committed Jun 15, 2017 1004 `````` void split_edges(Mesh * m, Node I, Edge a, Edge b, ElementArray & splitted_a, ElementArray & splitted_b, std::vector & transfer, bool print) `````` Kirill Terekhov committed Jul 18, 2016 1005 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1006 1007 `````` split_edge(m, I, a, splitted_a, transfer, print); split_edge(m, I, b, splitted_b, transfer, print); `````` Kirill Terekhov committed Jul 18, 2016 1008 `````` } `````` Kirill Terekhov committed Jul 28, 2016 1009 `````` `````` Kirill Terekhov committed Jun 21, 2017 1010 1011 `````` void intersect_naive(Mesh * m, ElementArray & segments, ElementArray & nodes, std::vector & transfer, Tag pnt, bool print) `````` Kirill Terekhov committed Jul 18, 2016 1012 `````` { `````` Kirill Terekhov committed Jul 28, 2016 1013 `````` //Tag pnt = m->CreateTag("PROJ_PNT"+m->GetLocalProcessorRank(),DATA_REAL,NODE,NODE,2); `````` Kirill Terekhov committed Aug 12, 2016 1014 `````` std::set intersections; `````` Kirill Terekhov committed Jun 15, 2017 1015 `````` std::vector initials(segments.size() * 2); `````` Kirill Terekhov committed Jul 28, 2016 1016 `````` MarkerType initial = m->CreatePrivateMarker(); `````` Kirill Terekhov committed Jul 18, 2016 1017 1018 `````` for (int k = 0; k < (int)segments.size(); ++k) { `````` Kirill Terekhov committed Jun 15, 2017 1019 1020 `````` initials[k * 2 + 0] = segments[k]->getBeg()->GetHandle(); initials[k * 2 + 1] = segments[k]->getEnd()->GetHandle(); `````` Kirill Terekhov committed Jul 28, 2016 1021 1022 `````` segments[k]->getBeg()->SetPrivateMarker(initial); segments[k]->getEnd()->SetPrivateMarker(initial); `````` Kirill Terekhov committed Jun 15, 2017 1023 1024 `````` intersections.insert(make_point(segments[k]->getBeg(), pnt)); intersections.insert(make_point(segments[k]->getEnd(), pnt)); `````` Kirill Terekhov committed Jul 18, 2016 1025 `````` } `````` Kirill Terekhov committed Jun 15, 2017 1026 `````` for (int i = 0; i < (int)segments.size(); ++i) `````` Kirill Terekhov committed Jul 18, 2016 1027 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1028 `````` for (int j = i + 1; j < (int)segments.size(); ++j) `````` Kirill Terekhov committed Jul 18, 2016 1029 `````` { `````` Kirill Terekhov committed Jun 21, 2017 1030 `````` std::pair I = intersect_segments(m, segments[i], segments[j], intersections, pnt, print); `````` Kirill Terekhov committed Jun 15, 2017 1031 `````` if (I.first) `````` Kirill Terekhov committed Jul 18, 2016 1032 1033 `````` { ElementArray splitted_a, splitted_b; `````` Kirill Terekhov committed Jun 21, 2017 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]; `````` Kirill Terekhov committed Jul 18, 2016 1037 1038 1039 1040 1041 `````` segments.push_back(splitted_a[1]); segments.push_back(splitted_b[1]); } } } `````` Kirill Terekhov committed Nov 03, 2016 1042 `````` //nodes.clear(); `````` Kirill Terekhov committed Jun 15, 2017 1043 `````` for (std::set::iterator it = intersections.begin(); it != intersections.end(); ++it) `````` Kirill Terekhov committed Jul 18, 2016 1044 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1045 `````` if (!m->GetPrivateMarker(it->node, initial)) `````` Kirill Terekhov committed Aug 12, 2016 1046 `````` nodes.push_back(it->node); `````` Kirill Terekhov committed Jul 18, 2016 1047 `````` } `````` Kirill Terekhov committed Jun 15, 2017 1048 1049 `````` for (int k = 0; k < (int)initials.size(); ++k) m->RemPrivateMarker(initials[k], initial); `````` Kirill Terekhov committed Jul 28, 2016 1050 `````` m->ReleasePrivateMarker(initial); `````` Kirill Terekhov committed Jul 18, 2016 1051 1052 `````` } `````` Kirill Terekhov committed Aug 12, 2016 1053 `````` void block_number_union(Element n, ElementArray & adj, Tag block, Tag write) `````` Kirill Terekhov committed Jul 08, 2016 1054 `````` { `````` Kirill Terekhov committed Aug 12, 2016 1055 `````` Storage::integer_array bn = n->IntegerArray(write); `````` Kirill Terekhov committed Jun 15, 2017 1056 `````` if (!adj.empty()) `````` Kirill Terekhov committed Jul 08, 2016 1057 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1058 1059 1060 `````` dynarray uni, tmp; //std::vector uni, tmp; for (ElementArray::size_type k = 0; k < adj.size(); ++k) `````` Kirill Terekhov committed Jul 08, 2016 1061 `````` { `````` Kirill Terekhov committed Aug 12, 2016 1062 `````` Storage::integer_array be = adj[k]->IntegerArray(block); `````` Kirill Terekhov committed Jun 15, 2017 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 committed Aug 12, 2016 1065 `````` uni.swap(tmp); `````` Kirill Terekhov committed Jul 08, 2016 1066 `````` } `````` Kirill Terekhov committed Jun 15, 2017 1067 `````` bn.replace(bn.begin(), bn.end(), uni.begin(), uni.end()); `````` Kirill Terekhov committed Jul 08, 2016 1068 `````` } `````` Kirill Terekhov committed Aug 12, 2016 1069 `````` else bn.clear(); `````` Kirill Terekhov committed Jul 08, 2016 1070 `````` } `````` Kirill Terekhov committed Jun 15, 2017 1071 `````` `````` Kirill Terekhov committed Aug 12, 2016 1072 1073 1074 `````` void block_number_union_merge(Element n, ElementArray & adj, Tag block, Tag write) { Storage::integer_array bn = n->IntegerArray(write); `````` Kirill Terekhov committed Jun 15, 2017 1075 `````` if (!adj.empty()) `````` Kirill Terekhov committed Jul 08, 2016 1076 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1077 1078 1079 `````` dynarray uni(bn.begin(), bn.end()), tmp(uni.size()); //std::vector uni(bn.begin(), bn.end()), tmp(uni.size()); for (ElementArray::size_type k = 0; k < adj.size(); ++k) `````` Kirill Terekhov committed Jul 08, 2016 1080 `````` { `````` Kirill Terekhov committed Aug 12, 2016 1081 `````` Storage::integer_array be = adj[k]->IntegerArray(block); `````` Kirill Terekhov committed Jun 15, 2017 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 committed Aug 12, 2016 1084 `````` uni.swap(tmp); `````` Kirill Terekhov committed Jul 08, 2016 1085 `````` } `````` Kirill Terekhov committed Jun 15, 2017 1086 `````` bn.replace(bn.begin(), bn.end(), uni.begin(), uni.end()); `````` Kirill Terekhov committed Jul 08, 2016 1087 `````` } `````` Kirill Terekhov committed Jul 11, 2016 1088 `````` } `````` Kirill Terekhov committed Aug 12, 2016 1089 `````` `````` Kirill Terekhov committed Jun 15, 2017 1090 `````` void block_number_intersection(Element n, const ElementArray & adj, Tag block, Tag write) `````` Kirill Terekhov committed Nov 04, 2016 1091 1092 `````` { Storage::integer_array bn = n->IntegerArray(write); `````` Kirill Terekhov committed Jun 15, 2017 1093 `````` if (!adj.empty()) `````` Kirill Terekhov committed Nov 04, 2016 1094 1095 `````` { Storage::integer_array be = adj[0]->IntegerArray(block); `````` Kirill Terekhov committed Jun 15, 2017 1096 1097 1098 `````` dynarray inter(be.begin(), be.end()), tmp(inter.size()); //std::vector inter(be.begin(), be.end()), tmp(inter.size()); for (ElementArray::size_type k = 1; k < adj.size(); ++k) `````` Kirill Terekhov committed Nov 04, 2016 1099 1100 `````` { be = adj[k]->IntegerArray(block); `````` Kirill Terekhov committed Jun 15, 2017 1101 `````` tmp.resize(std::set_intersection(inter.begin(), inter.end(), be.begin(), be.end(), tmp.begin()) - tmp.begin()); `````` Kirill Terekhov committed Nov 04, 2016 1102 1103 `````` inter.swap(tmp); } `````` Kirill Terekhov committed Jun 15, 2017 1104 `````` bn.replace(bn.begin(), bn.end(), inter.begin(), inter.end()); `````` Kirill Terekhov committed Nov 04, 2016 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()); } `````` Kirill Terekhov committed Jun 15, 2017 1111 1112 `````` struct mult_rec `````` Kirill Terekhov committed Jul 11, 2016 1113 `````` { `````` Kirill Terekhov committed Jun 15, 2017 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; }; `````` Kirill Terekhov committed Jul 08, 2016 1137 `````` `````` Kirill Terekhov committed Jun 21, 2017 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_recs; struct well_rec { welspecs_rec wspec; std::map wpi; std::map compdat; }; //typedef std::map compdat_wells; //typedef std::map welspecs_wells; typedef std::map wells_data; struct schedule { char tagnum; double value; HandleType element; }; typedef std::vector schedules; `````` Kirill Terekhov committed Jul 08, 2016 1202 1203 `````` void Mesh::LoadECL(std::string File) { `````` Kirill Terekhov committed Mar 26, 2019 1204 `````` bool parallel_read = true; `````` Kirill Terekhov committed May 28, 2017 1205 `````` char have_perm = 0; `````` Kirill Terekhov committed Jul 18, 2016 1206 `````` std::cout << std::scientific; `````` Kirill Terekhov committed Aug 14, 2017 1207 `````` int perform_splitting = 0; `````` Kirill Terekhov committed Aug 09, 2017 1208 1209 `````` bool project_perm = false; int split_degenerate = 0; `````` Kirill Terekhov committed Jun 15, 2017 1210 1211 1212 `````` bool check_topology = false; int verbosity = 0; for (INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); ++k) `````` Kirill Terekhov committed Aug 25, 2016 1213 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1214 `````` if (file_options[k].first == "ECL_SPLIT_GLUED") `````` Kirill Terekhov committed Aug 25, 2016 1215 `````` { `````` Kirill Terekhov committed Jan 18, 2019 1216 1217 1218 `````` if (file_options[k].second == "ALL") perform_splitting = 2; else if (file_options[k].second == "TRUE") `````` Kirill Terekhov committed Aug 14, 2017 1219 `````` perform_splitting = 1; `````` Kirill Terekhov committed Aug 25, 2016 1220 `````` else `````` Kirill Terekhov committed Aug 14, 2017 1221 `````` perform_splitting = 0; `````` Kirill Terekhov committed Aug 25, 2016 1222 `````` } `````` Kirill Terekhov committed Jun 15, 2017 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 committed Nov 04, 2016 1233 `````` { `````` Kirill Terekhov committed Jun 15, 2017 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 committed Nov 04, 2016 1242 1243 1244 1245 `````` project_perm = true; else project_perm = false; } `````` Kirill Terekhov committed Jun 15, 2017 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 committed Nov 23, 2015 1255 1256 1257 1258 `````` } std::vector old_nodes(NumberOfNodes()); { unsigned qq = 0; `````` Kirill Terekhov committed Jun 15, 2017 1259 `````` for (Mesh::iteratorNode it = BeginNode(); it != EndNode(); ++it) `````` Kirill Terekhov committed Nov 23, 2015 1260 1261 `````` old_nodes[qq++] = *it; } `````` Kirill Terekhov committed Jun 15, 2017 1262 1263 `````` if (!old_nodes.empty()) std::sort(old_nodes.begin(), old_nodes.end(), CentroidComparator(this)); `````` Kirill Terekhov committed Nov 04, 2016 1264 `````` `````` Kirill Terekhov committed Jun 15, 2017 1265 1266 `````` FILE * f = fopen(File.c_str(), "r"); if (f == NULL) `````` Kirill Terekhov committed Nov 04, 2016 1267 1268 1269 1270 `````` { std::cout << __FILE__ << ":" << __LINE__ << " cannot open file " << File << std::endl; throw BadFileName; } `````` Kirill Terekhov committed Jun 15, 2017 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, int> > fs(1, std::make_pair(std::make_pair(f, File), 0)); `````` Kirill Terekhov committed Oct 13, 2017 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 committed Nov 23, 2015 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; `````` Kirill Terekhov committed Jun 15, 2017 1287 `````` Storage::integer dims[3], mapaxis[6] = { 0, 1, 0, 0, 1, 0 }; `````` Kirill Terekhov committed Nov 23, 2015 1288 `````` Storage::real inrad = 0; `````` Kirill Terekhov committed Jun 21, 2017 1289 `````` std::vector xyz, perm, poro, tops, zcorn, ntg, pressure, pbub, swat, soil, sgas, thconr; `````` Kirill Terekhov committed Jun 15, 2017 1290 1291 1292 1293 1294 1295 1296 1297 1298 `````` std::vector actnum, satnum, eqlnum, pvtnum, rocknum; std::vector multiply; mult_rec mult_cur; std::vector editnnc; editnnc_rec editnnc_cur; std::vector faults; fault_rec fault_cur; std::vector multflt; multflt_rec multflt_cur; `````` Kirill Terekhov committed Jun 21, 2017 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 tsteps(1,0.0); //skip first tstep std::pair compdat_cur; std::pair welspecs_cur; std::pair wconprodinje_cur; `````` Kirill Terekhov committed Jun 15, 2017 1307 1308 `````` bool multiply_trans = false; bool multiply_props = false; `````` Kirill Terekhov committed Jun 21, 2017 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; `````` Kirill Terekhov committed Jun 15, 2017 1314 `````` while (!fs.empty()) `````` Kirill Terekhov committed Nov 23, 2015 1315 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1316 `````` while (fgets(readline, 2048, fs.back().first.first) != NULL) `````` Kirill Terekhov committed Nov 23, 2015 1317 1318 1319 `````` { fs.back().second++; //line number { `````` Kirill Terekhov committed Jun 15, 2017 1320 `````` if (readline[strlen(readline) - 1] == '\n') readline[strlen(readline) - 1] = '\0'; `````` Kirill Terekhov committed Jun 21, 2017 1321 1322 `````` UnStar(readline,readlines); //remove n* strcpy(readline,readlines); //copy back `````` Kirill Terekhov committed Nov 23, 2015 1323 `````` text_end = static_cast(strlen(readline)); `````` Kirill Terekhov committed Jun 15, 2017 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 committed Nov 23, 2015 1328 1329 `````` p = readline + text_start; pend = readline + text_end + 1; `````` Kirill Terekhov committed Oct 13, 2017 1330 1331 `````` strcpy(pupper,p); //upper case for keyword comparison for (char * q = pupper; q < pupper + (pend-p); q++) *q = toupper(*q); `````` Kirill Terekhov committed Nov 23, 2015 1332 `````` } `````` Kirill Terekhov committed Jun 15, 2017 1333 1334 1335 `````` if (p[0] == '-' && p[1] == '-') continue; //skip comment if (waitlines) { waitlines--; continue; } //skip meaningful lines switch (state) `````` Kirill Terekhov committed Nov 23, 2015 1336 1337 `````` { case ECL_NONE: `````` Kirill Terekhov committed Oct 13, 2017 1338 `````` if (!ECLSTRCMP(pupper, "END")) //end of data - don't read beyond `````` Kirill Terekhov committed Nov 23, 2015 1339 1340 1341 `````` { goto ecl_exit_loop; } `````` Kirill Terekhov committed Oct 13, 2017 1342 1343 `````` else if (!ECLSTRCMP(pupper, "INCLUDE")) state = ECL_INCLUDE; else if (!ECLSTRCMP(pupper, "BOX")) `````` Kirill Terekhov committed Nov 23, 2015 1344 `````` { ``````