mesh_ecl_file.cpp 200 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 ``````#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 `````` Kirill Terekhov committed Sep 28, 2020 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 ``````#define ECL_PERMXY 11 #define ECL_PERMXZ 12 #define ECL_PERMYZ 13 #define ECL_PORO 14 #define ECL_MAPAXIS 15 #define ECL_INRAD 16 #define ECL_COORDS 17 #define ECL_ZCORN 18 #define ECL_ACTNUM 19 #define ECL_SATNUM 20 #define ECL_THCONR 21 #define ECL_MULTIPLY 22 #define ECL_MULTIPLY_MUL 23 #define ECL_MULTIPLY_BLK 24 #define ECL_EDITNNC 25 #define ECL_EDITNNC_BLK 26 #define ECL_EDITNNC_MUL 27 #define ECL_NTG 28 #define ECL_FAULTS 29 #define ECL_FAULTS_BLK 30 #define ECL_FAULTS_DIR 31 #define ECL_MULTFLT 32 #define ECL_MULTFLT_MUL 33 #define ECL_PRESSURE 34 #define ECL_SWAT 35 #define ECL_SOIL 36 #define ECL_SGAS 37 #define ECL_PBUB 38 #define ECL_EQLNUM 39 #define ECL_ROCKNUM 40 #define ECL_PVTNUM 41 #define ECL_COMPDAT 42 #define ECL_COMPDAT_BLK 43 #define ECL_COMPDAT_OPEN 44 #define ECL_COMPDAT_SATNUM 45 #define ECL_COMPDAT_TRANS 46 #define ECL_COMPDAT_BORE 47 #define ECL_COMPDAT_PERM 48 #define ECL_COMPDAT_SKIN 49 #define ECL_COMPDAT_DFAC 50 #define ECL_COMPDAT_DIR 51 #define ECL_COMPDAT_RAD 52 #define ECL_DATES 53 #define ECL_DATES_MON 54 #define ECL_DATES_YEAR 55 #define ECL_DATES_HRS 56 #define ECL_WELSPECS 57 #define ECL_WELSPECS_GROUP 58 #define ECL_WELSPECS_I 59 #define ECL_WELSPECS_J 60 #define ECL_WELSPECS_Z 61 #define ECL_WELSPECS_PHASE 62 #define ECL_TSTEP 63 #define ECL_WCONPRODINJE 64 #define ECL_WCONPRODINJE_TYPE 65 #define ECL_WCONPRODINJE_OPEN 66 #define ECL_WCONPRODINJE_CTRL 67 #define ECL_WCONPRODINJE_RATES 68 #define ECL_WCONPRODINJE_BHP 69 `````` Kirill Terekhov committed Jun 15, 2017 168 `````` `````` Kirill Terekhov committed Nov 23, 2015 169 170 171 172 173 174 175 176 177 178 179 180 `````` #define ECL_GTYPE_NONE 0 #define ECL_GTYPE_TOPS 1 #define ECL_GTYPE_ZCORN 2 #define ECL_GTYPE_RADIAL 3 #define ECL_GTYPE_CARTESIAN 4 #define ECL_VAR_NONE 0 #define ECL_VAR_REAL 1 #define ECL_VAR_INT 2 `````` Kirill Terekhov committed Jul 08, 2016 181 ``````#define ECL_IJK_DATA(i,j,k) (i + ((j)+(k)*dims[1])*dims[0]) `````` Kirill Terekhov committed Jul 29, 2016 182 183 ``````#define ECL_IJK_COORDS(i,j,l,coord) ((((i) + (j)*(dims[0]+1))*2+(l))*3+(coord)) #define ECL_IJK_ZCORN(i,j,k,l) ((k)*dims[1]*dims[0]*8 + (1-(l)/4)*dims[1]*dims[0]*4 + (j)*dims[0]*4 + ((l)/2%2)*dims[0]*2 + (i)*2 + (l)%2) `````` Kirill Terekhov committed Jul 08, 2016 184 185 186 187 188 `````` //line sweep events #define SEG_START 0 #define SEG_END 1 `````` Kirill Terekhov committed May 28, 2017 189 190 191 ``````#define HAVE_PERM_X 1 #define HAVE_PERM_Y 2 #define HAVE_PERM_Z 4 `````` Kirill Terekhov committed Sep 28, 2020 192 193 194 ``````#define HAVE_PERM_XY 8 #define HAVE_PERM_XZ 16 #define HAVE_PERM_YZ 32 `````` Kirill Terekhov committed Nov 23, 2015 195 `````` `````` Kirill Terekhov committed Jun 15, 2017 196 `````` `````` Kirill Terekhov committed Jun 04, 2017 197 `````` `````` Kirill Terekhov committed Nov 23, 2015 198 199 ``````namespace INMOST { `````` Kirill Terekhov committed Jun 15, 2017 200 `````` `````` Kirill Terekhov committed Jun 21, 2017 201 `````` `````` Kirill Terekhov committed Feb 18, 2021 202 `````` void GetDXYZ(Element e, INMOST_DATA_REAL_TYPE & dx, INMOST_DATA_REAL_TYPE & dy, INMOST_DATA_REAL_TYPE & dz) `````` Kirill Terekhov committed Jun 21, 2017 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 `````` { Storage::real maxmin[6]; maxmin[0] = -1e20; maxmin[1] = 1e20; maxmin[2] = -1e20; maxmin[3] = 1e20; maxmin[4] = -1e20; maxmin[5] = 1e20; ElementArray 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 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 `````` void Faults(Mesh * m, MarkerType frac_markers) { //m->BeginModification(); const Storage::real mult = 0.95; Tag indexes = m->CreateTag("FRAC_TEMP_INDEXES", DATA_INTEGER, CELL, NONE, 2); //associate a node or a pair of nodes to each non-fracture edge that ends with a fracture node Tag cell2node = m->CreateTag("CELL2NODE", DATA_REFERENCE, CELL, NONE); //new faces to reconnect Tag connface = m->CreateTag("CONNFACE", DATA_REFERENCE, CELL, NONE); //Unmark any boundary faces for (Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) { if (!it->FrontCell().isValid() && it->GetMarker(frac_markers)) it->RemMarker(frac_markers); } //For each node create several new ones that will be used in open fracture for (Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) { ElementArray n_faces = it->getFaces(); //retrive all faces joining at the node INMOST_DATA_ENUM_TYPE num_frac_faces = 0; `````` Kirill Terekhov committed Feb 18, 2021 302 `````` for (size_t k = 0; k < n_faces.size(); ++k) if (n_faces[k].GetMarker(frac_markers)) num_frac_faces++; `````` Kirill Terekhov committed Jun 15, 2017 303 304 305 `````` if (num_frac_faces) { ElementArray n_cells = it->getCells(); //get all cells of the node `````` Kirill Terekhov committed Feb 18, 2021 306 `````` for (size_t k = 0; k < n_cells.size(); ++k) //assign local enumeration to the cells `````` Kirill Terekhov committed Feb 21, 2021 307 `````` n_cells[k].IntegerDF(indexes) = static_cast(k); `````` Kirill Terekhov committed Feb 18, 2021 308 `````` for (size_t k = 0; k < n_faces.size(); ++k) //stich together node's numbers corresponding to cells if no fracture separates them `````` Kirill Terekhov committed Jun 15, 2017 309 310 311 312 313 314 315 `````` { if (!n_faces[k].GetMarker(frac_markers) && n_faces[k].FrontCell().isValid()) { int bi = n_faces[k].BackCell()->IntegerDF(indexes); int fi = n_faces[k].FrontCell()->IntegerDF(indexes); if (bi != fi) { `````` Kirill Terekhov committed Feb 18, 2021 316 `````` for (size_t q = 0; q < n_cells.size(); q++) if (n_cells[q].IntegerDF(indexes) == fi) `````` Kirill Terekhov committed Jun 15, 2017 317 318 319 320 321 `````` n_cells[q].IntegerDF(indexes) = bi; } } } dynarray nums(n_cells.size()); //store all numbers `````` Kirill Terekhov committed Feb 18, 2021 322 `````` for (size_t k = 0; k < n_cells.size(); ++k) nums[k] = n_cells[k].IntegerDF(indexes); `````` Kirill Terekhov committed Jun 15, 2017 323 324 325 326 `````` std::sort(nums.begin(), nums.end()); nums.resize(std::unique(nums.begin(), nums.end()) - nums.begin()); if (nums.size() > 1) //at least two distinctive nodes { `````` Kirill Terekhov committed Feb 18, 2021 327 `````` for (size_t k = 0; k < nums.size(); k++) `````` Kirill Terekhov committed Jun 15, 2017 328 329 330 331 `````` { Node image; Storage::real xyz[3] = { 0, 0, 0 }, cntc[3] = { 0, 0, 0 }; int n_node_cells = 0; `````` Kirill Terekhov committed Feb 18, 2021 332 `````` for (size_t q = 0; q < n_cells.size(); q++) `````` Kirill Terekhov committed Jun 15, 2017 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 `````` { if (n_cells[q].IntegerDF(indexes) == nums[k]) { n_cells[q].Centroid(cntc); xyz[0] += cntc[0]; xyz[1] += cntc[1]; xyz[2] += cntc[2]; n_node_cells++; } } xyz[0] /= static_cast(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); `````` Kirill Terekhov committed Feb 18, 2021 351 `````` for (size_t q = 0; q < n_cells.size(); q++) `````` Kirill Terekhov committed Jun 15, 2017 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 `````` { if (n_cells[q].IntegerDF(indexes) == nums[k]) { Storage::reference_array fnodes = n_cells[q].ReferenceArray(cell2node); fnodes.push_back(*it); //TODO fnodes.push_back(image); } } } it->SetMarker(frac_markers); //mark node as fracture node //std::cout << "Number of nodes: " << m->NumberOfNodes() << std::endl; } } } m->DeleteTag(indexes); //Mark edges that will be incorporated into fracture cells for (Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it) { if (it->nbAdjElements(FACE, frac_markers) > 2) it->SetMarker(frac_markers); } //create new matrix-matrix faces, adjacent to fractures //for all non-fracture faces that have any fracture edge create new face for (Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if (!it->GetMarker(frac_markers) && it->nbAdjElements(NODE, frac_markers)) { ElementArray 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 `````` Kirill Terekhov committed Feb 18, 2021 432 `````` for (size_t q = 0; q < nodes.size(); ++q) if (nodes[q].GetMarker(frac_markers)) `````` Kirill Terekhov committed Jun 15, 2017 433 434 435 436 437 438 439 440 441 `````` { Node n; for (int k = 0; k < (int)images.size() && !n.isValid(); k += 2) if (images[k] == nodes[q]) n = images[k + 1]->getAsNode(); assert(n.isValid()); nodesb[q] = n; } else nodesb[q] = nodes[q]; //create edges `````` Kirill Terekhov committed Feb 18, 2021 442 `````` for (size_t k = 0; k < edges.size(); k++) `````` Kirill Terekhov committed Jun 15, 2017 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 `````` { if (nodes[k].GetMarker(frac_markers) || nodes[(k + 1) % N].GetMarker(frac_markers)) { ElementArray 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 `````` Kirill Terekhov committed Feb 18, 2021 464 `````` for (size_t q = 0; q < nodes.size(); ++q) if (nodes[q].GetMarker(frac_markers)) `````` Kirill Terekhov committed Jun 15, 2017 465 466 467 468 469 470 471 472 473 `````` { Node n; for (int k = 0; k < (int)images.size() && !n.isValid(); k += 2) if (images[k] == nodes[q]) n = images[k + 1]->getAsNode(); assert(n.isValid()); nodesf[q] = n; } else nodesf[q] = nodes[q]; //create edges `````` Kirill Terekhov committed Feb 18, 2021 474 `````` for (size_t k = 0; k < edges.size(); k++) `````` Kirill Terekhov committed Jun 15, 2017 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 `````` { if (nodes[k].GetMarker(frac_markers) || nodes[(k + 1) % N].GetMarker(frac_markers)) { ElementArray 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 531 `````` //replace n* witn n stars `````` Kirill Terekhov committed Feb 18, 2021 532 `````` void UnStar(const char * input, char * output) `````` Kirill Terekhov committed Jun 21, 2017 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 `````` { size_t q = 1, j; output[0] = input[0]; //we skip first char if (input[0] != '\0') { for (size_t k = 1; k < strlen(input) - 1; ++k) { //current character is star, next char is space or end of line and prev char is not space if (input[k] == '*' && (isspace(input[k + 1]) || input[k + 1] == '/') && !isspace(input[k - 1])) { j = k; //roll back to previous space while (!isspace(input[k]) && k > 0) { k--; q--; } //advance from space if (isspace(input[k])) { k++; q++; } int count = atoi(input+k); for (int l = 0; l < count; ++l) { output[q++] = '*'; output[q++] = ' '; } k = j; } else output[q++] = input[k]; } output[q++] = input[strlen(input) - 1]; //we skip last char output[q++] = '\0'; // end of line } } `````` Kirill Terekhov committed Feb 18, 2021 571 `````` std::string STR8(const char * input) `````` Kirill Terekhov committed Jun 15, 2017 572 573 574 575 576 577 578 579 `````` { std::string ret(8, ' '); int kend = (int)strlen(input); if (kend > 8) kend = 8; for (int k = 0; k < kend; k++) ret[k] = input[k]; return ret; } `````` Kirill Terekhov committed Feb 18, 2021 580 `````` std::string UnQoute(std::string input) `````` Kirill Terekhov committed Jun 15, 2017 581 582 583 584 585 586 587 `````` { if (input[0] == '\'' && input[input.size() - 1] == '\'') return input.substr(1, input.size() - 2); else if (input[0] == '"' && input[input.size() - 1] == '"') return input.substr(1, input.size() - 2); return input; } `````` Kirill Terekhov committed Feb 18, 2021 588 `````` std::string ToUpper(std::string input) `````` Kirill Terekhov committed Jun 15, 2017 589 590 591 592 593 594 `````` { for (size_t k = 0; k < input.size(); ++k) input[k] = toupper(input[k]); return input; } `````` Kirill Terekhov committed Feb 18, 2021 595 `````` int GetDir(std::string dir) `````` Kirill Terekhov committed Jun 15, 2017 596 597 598 `````` { if (dir == "X-" || dir == "I-") return 0; `````` Kirill Terekhov committed Aug 05, 2017 599 `````` else if (dir == "X" || dir == "I" || dir == "X+" || dir == "I+" || dir == "X " || dir == "I ") `````` Kirill Terekhov committed Jun 15, 2017 600 601 602 `````` return 1; else if (dir == "Y-" || dir == "J-") return 2; `````` Kirill Terekhov committed Aug 05, 2017 603 `````` else if (dir == "Y" || dir == "J" || dir == "Y+" || dir == "J+" || dir == "Y " || dir == "J ") `````` Kirill Terekhov committed Jun 15, 2017 604 605 606 `````` return 3; else if (dir == "Z-" || dir == "K-") return 4; `````` Kirill Terekhov committed Aug 05, 2017 607 `````` else if (dir == "Z" || dir == "K" || dir == "Z+" || dir == "K+" || dir == "Z " || dir == "K ") `````` Kirill Terekhov committed Jun 15, 2017 608 609 610 `````` return 5; return -1; } `````` Kirill Terekhov committed Feb 18, 2021 611 `````` std::string GetFolder(std::string file) `````` Kirill Terekhov committed Jul 27, 2016 612 613 `````` { size_t found = file.find_last_of("/\\"); `````` Kirill Terekhov committed Jun 15, 2017 614 `````` if (found == std::string::npos) `````` Kirill Terekhov committed Jul 27, 2016 615 `````` return ""; `````` Kirill Terekhov committed Jun 15, 2017 616 617 618 `````` else return file.substr(0, found); } `````` Kirill Terekhov committed Feb 18, 2021 619 `````` ElementArray OrderEdges(const ElementArray & input, MarkerType mrk) `````` Kirill Terekhov committed Jun 15, 2017 620 621 622 623 `````` { Mesh * m = const_cast(input.GetMeshLink()); ElementArray ret(m); ret.reserve(input.size()); `````` Kirill Terekhov committed Feb 19, 2021 624 `````` HandleType last = InvalidHandle(); `````` Kirill Terekhov committed Jun 15, 2017 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 `````` ret.push_back(input[0]); HandleType add = InvalidHandle(); const Element::adj_type & lc0 = m->LowConn(input.at(0)); for (int l = 1; l < (int)input.size(); ++l) { const Element::adj_type & lcl = m->LowConn(input.at(l)); if (lcl[0] == lc0[0] || lcl[0] == lc0[1]) { add = input.at(l); last = lcl[1]; break; } else if (lcl[1] == lc0[0] || lcl[1] == lc0[1]) { add = input.at(l); last = lcl[0]; break; } } if (add == InvalidHandle()) std::cout << __FILE__ << ":" << __LINE__ << " bad edges, number of edges " << input.size() << std::endl; else { ret.push_back(add); m->SetPrivateMarker(add, mrk); } while (ret.size() != input.size()) { add = InvalidHandle(); for (int l = 1; l < (int)input.size(); ++l) { if (!m->GetPrivateMarker(input.at(l), mrk)) { const Element::adj_type & lcl = m->LowConn(input.at(l)); if (lcl[0] == last) { add = input.at(l); last = lcl[1]; break; } else if (lcl[1] == last) { add = input.at(l); last = lcl[0]; break; } } } if (add == InvalidHandle()) std::cout << __FILE__ << ":" << __LINE__ << " bad edges, number of edges " << input.size() << std::endl; else { ret.push_back(add); m->SetPrivateMarker(add, mrk); } } ret.RemPrivateMarker(mrk); return ret; `````` Kirill Terekhov committed Jul 27, 2016 684 `````` } `````` Kirill Terekhov committed Jun 15, 2017 685 `````` `````` Kirill Terekhov committed Jun 21, 2017 686 `````` `````` Kirill Terekhov committed Jul 08, 2016 687 688 689 `````` //2d point with comparison operator for line sweep algorithm struct Point { `````` Kirill Terekhov committed Aug 12, 2016 690 `````` HandleType node; `````` Kirill Terekhov committed Nov 04, 2016 691 `````` Storage::real x, y; `````` Kirill Terekhov committed Aug 12, 2016 692 693 694 695 `````` Point & operator = (Point const & b) { node = b.node; x = b.x; y = b.y; return *this; } Point(const Point & b) : node(b.node), x(b.x), y(b.y) {} Point(HandleType _node, Storage::real _x, Storage::real _y) : node(_node), x(_x), y(_y) {} Point(HandleType _node, Storage::real v[2]) : node(_node), x(v[0]), y(v[1]) {} `````` Kirill Terekhov committed Jul 08, 2016 696 697 `````` bool operator <(const Point & b) const { `````` Kirill Terekhov committed Aug 12, 2016 698 699 700 `````` if (y < b.y - ECL_POINT_EPS) return true; else if (y > b.y + ECL_POINT_EPS) return false; else if (x < b.x - ECL_POINT_EPS) return true; `````` Kirill Terekhov committed Jul 08, 2016 701 702 703 704 `````` else return false; } bool operator ==(const Point & b) const { `````` Kirill Terekhov committed Aug 12, 2016 705 `````` return std::abs(y - b.y) < ECL_POINT_EPS && std::abs(x - b.x) < ECL_POINT_EPS; `````` Kirill Terekhov committed Jul 08, 2016 706 707 708 `````` } bool operator !=(const Point & b) const { `````` Kirill Terekhov committed Aug 12, 2016 709 `````` return std::abs(y - b.y) > ECL_POINT_EPS || std::abs(x - b.x) > ECL_POINT_EPS; `````` Kirill Terekhov committed Jul 29, 2016 710 711 `````` } }; `````` Kirill Terekhov committed Jun 15, 2017 712 `````` `````` Kirill Terekhov committed Aug 12, 2016 713 `````` // Used to sort arrays by indices `````` Kirill Terekhov committed Jul 18, 2016 714 715 716 717 718 719 `````` class index_comparator { Storage::integer_array & data; public: index_comparator(Storage::integer_array & data) : data(data) {} index_comparator(const index_comparator & b) : data(b.data) {} `````` Kirill Terekhov committed Jun 15, 2017 720 721 `````` index_comparator & operator =(index_comparator const & b) { data = b.data; return *this; } bool operator ()(int a, int b) const { return data[a] < data[b]; } `````` Kirill Terekhov committed Jul 18, 2016 722 `````` }; `````` Kirill Terekhov committed Aug 12, 2016 723 724 725 726 727 728 729 `````` // Used to sort arrays by depth class depth_comparator { const Storage::real * data; public: depth_comparator(const Storage::real * data) : data(data) {} depth_comparator(const depth_comparator & b) : data(b.data) {} `````` Kirill Terekhov committed Jun 15, 2017 730 731 `````` depth_comparator & operator =(depth_comparator const & b) { data = b.data; return *this; } bool operator ()(int a, int b) const { return data[a] < data[b]; } `````` Kirill Terekhov committed Aug 12, 2016 732 `````` }; `````` Kirill Terekhov committed Jun 15, 2017 733 `````` `````` Kirill Terekhov committed Jul 18, 2016 734 735 736 737 `````` template int count_duplicates(ElementArray & array) { Mesh * m = array.GetMeshLink(); `````` Kirill Terekhov committed Jul 28, 2016 738 `````` MarkerType mrk = m->CreatePrivateMarker(); `````` Kirill Terekhov committed Jul 18, 2016 739 `````` int dups = 0; `````` Kirill Terekhov committed Feb 18, 2021 740 `````` for (typename ElementArray::size_type k = 0; k < array.size(); ++k) `````` Kirill Terekhov committed Jul 18, 2016 741 `````` { `````` Kirill Terekhov committed Jun 15, 2017 742 `````` if (array[k].GetPrivateMarker(mrk)) `````` Kirill Terekhov committed Jul 18, 2016 743 `````` dups++; `````` Kirill Terekhov committed Jul 28, 2016 744 `````` array[k].SetPrivateMarker(mrk); `````` Kirill Terekhov committed Jul 18, 2016 745 `````` } `````` Kirill Terekhov committed Jul 28, 2016 746 747 `````` array.RemPrivateMarker(mrk); m->ReleasePrivateMarker(mrk); `````` Kirill Terekhov committed Jul 18, 2016 748 749 `````` return dups; } `````` Kirill Terekhov committed Jun 15, 2017 750 `````` `````` Kirill Terekhov committed Jul 18, 2016 751 752 753 754 `````` template void make_unique(ElementArray & array) { Mesh * m = array.GetMeshLink(); `````` Kirill Terekhov committed Jul 28, 2016 755 `````` MarkerType mrk = m->CreatePrivateMarker(); `````` Kirill Terekhov committed Feb 18, 2021 756 `````` typename ElementArray::size_type i = 0, j = 0; `````` Kirill Terekhov committed Jun 15, 2017 757 `````` while (j < array.size()) `````` Kirill Terekhov committed Jul 18, 2016 758 `````` { `````` Kirill Terekhov committed Jun 15, 2017 759 `````` if (array[j].GetPrivateMarker(mrk)) `````` Kirill Terekhov committed Jul 18, 2016 760 761 762 `````` ++j; else { `````` Kirill Terekhov committed Jul 28, 2016 763 `````` array[j].SetPrivateMarker(mrk); `````` Kirill Terekhov committed Jul 18, 2016 764 765 766 `````` array[i++] = array[j++]; } } `````` Kirill Terekhov committed Jul 28, 2016 767 `````` array.RemPrivateMarker(mrk); `````` Kirill Terekhov committed Jul 18, 2016 768 `````` array.resize(i); `````` Kirill Terekhov committed Jul 28, 2016 769 `````` m->ReleasePrivateMarker(mrk); `````` Kirill Terekhov committed Jul 18, 2016 770 `````` } `````` Kirill Terekhov committed Jun 15, 2017 771 `````` `````` Kirill Terekhov committed Aug 12, 2016 772 773 `````` Point make_point(Node n, Tag pnt) { `````` Kirill Terekhov committed Jun 15, 2017 774 `````` return Point(n->GetHandle(), n->RealArray(pnt).data()); `````` Kirill Terekhov committed Aug 12, 2016 775 `````` } `````` Kirill Terekhov committed Nov 04, 2016 776 `````` `````` Kirill Terekhov committed Jun 15, 2017 777 `````` std::pair intersect_pairs(const Point & pabeg, const Point & paend, const Point & pbbeg, const Point & pbend) `````` Kirill Terekhov committed Nov 04, 2016 778 779 `````` { const Storage::real eps = ECL_INTERSECT_EPS; `````` Kirill Terekhov committed Jun 15, 2017 780 `````` Point pfind(InvalidHandle(), 0, 0); `````` Kirill Terekhov committed Feb 19, 2021 781 `````` Storage::real t1 = -1, t2 = -1; `````` Kirill Terekhov committed Nov 04, 2016 782 783 `````` Storage::real div = (pabeg.x - paend.x)*(pbbeg.y - pbend.y) - (pabeg.y - paend.y)*(pbbeg.x - pbend.x); if (std::abs(div) < 1.0e-50) `````` Kirill Terekhov committed Jun 15, 2017 784 `````` return std::make_pair(-1, -1); `````` Kirill Terekhov committed Nov 04, 2016 785 786 787 `````` pfind.x = ((pabeg.x*paend.y - pabeg.y*paend.x)*(pbbeg.x - pbend.x) - (pabeg.x - paend.x)*(pbbeg.x*pbend.y - pbbeg.y*pbend.x)) / div; pfind.y = ((pabeg.x*paend.y - pabeg.y*paend.x)*(pbbeg.y - pbend.y) - (pabeg.y - paend.y)*(pbbeg.x*pbend.y - pbbeg.y*pbend.x)) / div; //optimization! - should belong unit cube `````` Kirill Terekhov committed Jun 15, 2017 788 `````` 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 789 790 791 `````` if (std::abs(paend.x - pabeg.x) > eps) { t1 = (pfind.x - pabeg.x) / (paend.x - pabeg.x); `````` Kirill Terekhov committed Jun 15, 2017 792 `````` if (t1 < eps || t1 > 1.0 - eps) { return std::make_pair(-1, -1); } `````` Kirill Terekhov committed Nov 04, 2016 793 794 795 796 `````` } if (std::abs(paend.y - pabeg.y) > eps) { t1 = (pfind.y - pabeg.y) / (paend.y - pabeg.y); `````` Kirill Terekhov committed Jun 15, 2017 797 `````` if (t1 < eps || t1 > 1.0 - eps) { return std::make_pair(-1, -1); } `````` Kirill Terekhov committed Nov 04, 2016 798 799 800 801 `````` } if (std::abs(pbend.x - pbbeg.x) > eps) { t2 = (pfind.x - pbbeg.x) / (pbend.x - pbbeg.x); `````` Kirill Terekhov committed Jun 15, 2017 802 `````` if (t2 < eps || t2 > 1.0 - eps) { return std::make_pair(-1, -1); } `````` Kirill Terekhov committed Nov 04, 2016 803 804 805 806 `````` } if (std::abs(pbend.y - pbbeg.y) > eps) { t2 = (pfind.y - pbbeg.y) / (pbend.y - pbbeg.y); `````` Kirill Terekhov committed Jun 15, 2017 807 `````` if (t2 < eps || t2 > 1.0 - eps) { return std::make_pair(-1, -1); } `````` Kirill Terekhov committed Nov 04, 2016 808 `````` } `````` Kirill Terekhov committed Jun 15, 2017 809 `````` return std::make_pair(t1, t2); `````` Kirill Terekhov committed Nov 04, 2016 810 `````` } `````` Kirill Terekhov committed Jun 15, 2017 811 `````` `````` Kirill Terekhov committed Jun 21, 2017 812 813 814 `````` //returns distance between line if shortest line is within segments, writes position of distance into last argument `````` Kirill Terekhov committed Feb 18, 2021 815 `````` static INMOST_DATA_REAL_TYPE SegmentDistance(const INMOST_DATA_REAL_TYPE v1[3], const INMOST_DATA_REAL_TYPE v2[3], const INMOST_DATA_REAL_TYPE v3[3], const INMOST_DATA_REAL_TYPE v4[3], INMOST_DATA_REAL_TYPE vout[3]) `````` Kirill Terekhov committed Jun 21, 2017 816 `````` { `````` Kirill Terekhov committed Feb 18, 2021 817 `````` INMOST_DATA_REAL_TYPE v13[3], v21[3], v43[3]; `````` Kirill Terekhov committed Jun 21, 2017 818 819 820 821 822 823 `````` for (int k = 0; k < 3; ++k) { v13[k] = v1[k] - v3[k]; v21[k] = v2[k] - v1[k]; v43[k] = v4[k] - v3[k]; } `````` Kirill Terekhov committed Feb 18, 2021 824 `````` INMOST_DATA_REAL_TYPE d1321 = 0, d2121 = 0, d4321 = 0, d1343 = 0, d4343 = 0; `````` Kirill Terekhov committed Jun 21, 2017 825 826 827 828 829 830 831 832 `````` for (int k = 0; k < 3; ++k) { d1321 += v13[k] * v21[k]; d2121 += v21[k] * v21[k]; d4321 += v43[k] * v21[k]; d1343 += v13[k] * v43[k]; d4343 += v43[k] * v43[k]; } `````` Kirill Terekhov committed Feb 18, 2021 833 `````` INMOST_DATA_REAL_TYPE mu1 = 0, mu2 = 0; `````` Kirill Terekhov committed Jun 21, 2017 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 `````` /* double parallel = 0; for (int k = 0; k < 3; ++k) parallel += v21[k] * v43[k]; parallel = fabs(parallel); parallel /= sqrt(d2121)*sqrt(d4343); if (fabs(parallel - 1.0) > 1.0e-6) */ if (fabs(d2121*d4343 - d4321*d4321) > 1.0e-6) { mu1 = (d1343*d4321 - d1321*d4343) / (d2121*d4343 - d4321*d4321); mu2 = (d1343 + mu1*d4321) / d4343; if (mu1 > 1) mu1 = 1; if (mu1 < 0) mu1 = 0; if (mu2 > 1) mu2 = 1; if (mu2 < 0) mu2 = 0; } else //parallel lines { //1d problem `````` Kirill Terekhov committed Feb 18, 2021 854 855 `````` INMOST_DATA_REAL_TYPE la = 0, ra = 0; INMOST_DATA_REAL_TYPE lb = 0, rb = 0; `````` Kirill Terekhov committed Jun 21, 2017 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 `````` for (int k = 0; k < 3; ++k) { la += v1[k] * v21[k]; ra += v2[k] * v21[k]; lb += v3[k] * v21[k]; rb += v4[k] * v21[k]; } bool invb = false; if (lb > rb) { std::swap(lb, rb); invb = true; } if (ra+1.0e-7 < lb) // (la,ra) is to the left of (lb,rb), no intersection { mu1 = 1; mu2 = 0; } else if (la > rb+1.0e-7) // (la,ra) is to the right of (lb,rb), no intersection { mu1 = 0; mu2 = 1; } else if (lb > la) // (lb,rb) intersects to the right of or contains (la,ra) { mu2 = 0; mu1 = (lb - la) / (ra - la); } else if (rb < ra) // (lb,rb) intersects to the left of or contains (la,ra) { mu2 = 1; mu1 = (rb - la) / (ra - la); } else // (lb,rb) contains (la,ra) { mu1 = 0; mu2 = (la - lb) / (rb - lb); } if (invb) mu2 = 1 - mu2; } `````` Kirill Terekhov committed Feb 18, 2021 897 `````` INMOST_DATA_REAL_TYPE vs1[3], vs2[3], h = 0; `````` Kirill Terekhov committed Jun 21, 2017 898 899 900 901 902 903 904 905 906 907 908 909 910 `````` for (int k = 0; k < 3; ++k) { vs1[k] = v1[k] + mu1*(v2[k] - v1[k]); vs2[k] = v3[k] + mu2*(v4[k] - v3[k]); vout[k] = (vs1[k] + vs2[k])*0.5; h += (vs2[k] - vs1[k])*(vs2[k] - vs1[k]); } return sqrt(h); } std::pair intersect_segments(Mesh * m, const Edge & a, const Edge & b, std::set & intersections, Tag pnt, bool print) `````` Kirill Terekhov committed Jul 08, 2016 911 `````` { `````` Kirill Terekhov committed Aug 12, 2016 912 `````` const Storage::real eps = ECL_INTERSECT_EPS; `````` Kirill Terekhov committed Jun 15, 2017 913 914 915 916 917 918 919 920 921 922 `````` Point pfind(InvalidHandle(), 0, 0); if (a->getBeg() == b->getBeg() || a->getBeg() == b->getEnd() || a->getEnd() == b->getEnd() || a->getEnd() == b->getBeg()) return std::make_pair(false, pfind); Point pabeg = make_point(a->getBeg(), pnt); Point paend = make_point(a->getEnd(), pnt); Point pbbeg = make_point(b->getBeg(), pnt); Point pbend = make_point(b->getEnd(), pnt); `````` Kirill Terekhov committed Nov 04, 2016 923 `````` 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 924 `````` Storage::real t1, t2; `````` Kirill Terekhov committed Aug 12, 2016 925 `````` if (std::abs(div) < 1.0e-50) `````` Kirill Terekhov committed Jul 18, 2016 926 927 `````` { if (print) std::cout << "divisor is zero" << std::endl; `````` Kirill Terekhov committed Jun 15, 2017 928 `````` return std::make_pair(false, pfind); `````` Kirill Terekhov committed Jul 18, 2016 929 930 931 `````` } pfind.x = ((pabeg.x*paend.y - pabeg.y*paend.x)*(pbbeg.x - pbend.x) - (pabeg.x - paend.x)*(pbbeg.x*pbend.y - pbbeg.y*pbend.x)) / div; pfind.y = ((pabeg.x*paend.y - pabeg.y*paend.x)*(pbbeg.y - pbend.y) - (pabeg.y - paend.y)*(pbbeg.x*pbend.y - pbbeg.y*pbend.x)) / div; `````` Kirill Terekhov committed Jul 28, 2016 932 `````` //optimization uses information that we stay in unit cube `````` Kirill Terekhov committed Jun 15, 2017 933 934 935 `````` if (pfind.x < 0 || pfind.x > 1 || pfind.y < 0 || pfind.y > 1) return std::make_pair(false, pfind); if (print) std::cout << "found (" << pfind.x << ", " << pfind.y << ") for edges " << a->GetHandle() << " and " << b->GetHandle() << std::endl; `````` Kirill Terekhov committed Jul 18, 2016 936 `````` //probably some of these tests are redundant `````` Kirill Terekhov committed Aug 12, 2016 937 `````` if (std::abs(paend.x - pabeg.x) > eps) `````` Kirill Terekhov committed Jul 18, 2016 938 939 `````` { t1 = (pfind.x - pabeg.x) / (paend.x - pabeg.x); `````` Kirill Terekhov committed Jun 15, 2017 940 `````` 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 941 `````` } `````` Kirill Terekhov committed Aug 12, 2016 942 `````` if (std::abs(paend.y - pabeg.y) > eps) `````` Kirill Terekhov committed Jul 18, 2016 943 944 `````` { t1 = (pfind.y - pabeg.y) / (paend.y - pabeg.y); `````` Kirill Terekhov committed Jun 15, 2017 945 `````` 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 946 `````` } `````` Kirill Terekhov committed Aug 12, 2016 947 `````` if (std::abs(pbend.x - pbbeg.x) > eps) `````` Kirill Terekhov committed Jul 18, 2016 948 949 `````` { t2 = (pfind.x - pbbeg.x) / (pbend.x - pbbeg.x); `````` Kirill Terekhov committed Jun 15, 2017 950 `````` 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 951 `````` } `````` Kirill Terekhov committed Aug 12, 2016 952 `````` if (std::abs(pbend.y - pbbeg.y) > eps) `````` Kirill Terekhov committed Jul 18, 2016 953 954 `````` { t2 = (pfind.y - pbbeg.y) / (pbend.y - pbbeg.y); `````` Kirill Terekhov committed Jun 15, 2017 955 `````` 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 956 `````` } `````` Kirill Terekhov committed Aug 12, 2016 957 `````` std::set::iterator search = intersections.find(pfind); `````` Kirill Terekhov committed Jul 18, 2016 958 `````` //check whether intersection already exists `````` Kirill Terekhov committed Jun 15, 2017 959 `````` if (search != intersections.end()) //there is a node! `````` Kirill Terekhov committed Aug 12, 2016 960 961 `````` { pfind = *search; `````` Kirill Terekhov committed Jun 15, 2017 962 963 964 965 966 `````` if (pfind.node == a->getBeg()->GetHandle() || pfind.node == b->getBeg()->GetHandle() || pfind.node == a->getEnd()->GetHandle() || pfind.node == b->getEnd()->GetHandle()) return std::make_pair(false, pfind); `````` Kirill Terekhov committed Aug 12, 2016 967 `````` } `````` Kirill Terekhov committed Jul 18, 2016 968 969 `````` else //no node, create one { `````` Kirill Terekhov committed Aug 12, 2016 970 `````` Storage::real find[3]; `````` Kirill Terekhov committed Jul 29, 2016 971 `````` //restore coordinate `````` Kirill Terekhov committed Jun 21, 2017 972 973 `````` SegmentDistance(a->getBeg().Coords().data(), a->getEnd().Coords().data(), b->getBeg().Coords().data(), b->getEnd().Coords().data(), find); `````` Kirill Terekhov committed Aug 12, 2016 974 975 `````` pfind.node = m->CreateNode(find)->GetHandle(); if (print) std::cout << "intersection accepted (" << find[0] << "," << find[1] << "," << find[2] << ") t1 " << t1 << " t2 " << t2 << " new node " << pfind.node << std::endl; `````` Kirill Terekhov committed Jun 15, 2017 976 `````` Storage::real_array _pfind = m->RealArray(pfind.node, pnt); `````` Kirill Terekhov committed Jul 28, 2016 977 978 `````` _pfind[0] = pfind.x; _pfind[1] = pfind.y; `````` Kirill Terekhov committed Aug 12, 2016 979 `````` intersections.insert(pfind); `````` Kirill Terekhov committed Jul 18, 2016 980 `````` } `````` Kirill Terekhov committed Jun 15, 2017 981 `````` return std::make_pair(true, pfind); `````` Kirill Terekhov committed Jul 18, 2016 982 `````` } `````` Kirill Terekhov committed Jul 08, 2016 983 `````` `````` Kirill Terekhov committed Nov 03, 2016 984 985 986 987 `````` 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 Feb 18, 2021 988 `````` for (size_t k = 0; k < transfer.size(); ++k) //tags `````` Kirill Terekhov committed Nov 03, 2016 989 `````` { `````` Kirill Terekhov committed Feb 18, 2021 990 `````` INMOST_DATA_ENUM_TYPE size = a->GetDataSize(transfer[k]); `````` Kirill Terekhov committed Feb 21, 2021 991 `````` copy[k].resize((size_t)transfer[k].GetBytesSize()*size); `````` Kirill Terekhov committed Jun 15, 2017 992 `````` if (!copy.empty()) a->GetData(transfer[k], 0, size, ©[k][0]); `````` Kirill Terekhov committed Nov 03, 2016 993 `````` } `````` Kirill Terekhov committed Jun 15, 2017 994 995 996 997 `````` 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 998 `````` //duplicate data `````` Kirill Terekhov committed Feb 18, 2021 999 `````` for (size_t k = 0; k < transfer.size(); ++k) `````` Kirill Terekhov committed Nov 03, 2016 1000 `````` { `````` Kirill Terekhov committed Feb 21, 2021 1001 `````` INMOST_DATA_ENUM_TYPE size = static_cast(copy[k].size() / transfer[k].GetBytesSize()); `````` Kirill Terekhov committed Jun 15, 2017 1002 `````` if (size) for (int l = 0; l < 2; ++l) //two parts `````` Kirill Terekhov committed Nov 03, 2016 1003 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1004 1005 `````` splitted_a[l].SetDataSize(transfer[k], size); splitted_a[l].SetData(transfer[k], 0, size, ©[k][0]); `````` Kirill Terekhov committed Nov 03, 2016 1006 1007 1008 `````` } } } `````` Kirill Terekhov committed Jul 18, 2016 1009 `````` `````` Kirill Terekhov committed Jun 15, 2017 1010 `````` 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 1011 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1012 1013 `````` split_edge(m, I, a, splitted_a, transfer, print); split_edge(m, I, b, splitted_b, transfer, print); `````` Kirill Terekhov committed Jul 18, 2016 1014 `````` } `````` Kirill Terekhov committed Jul 28, 2016 1015 `````` `````` Kirill Terekhov committed Jun 21, 2017 1016 1017 `````` void intersect_naive(Mesh * m, ElementArray & segments, ElementArray & nodes, std::vector & transfer, Tag pnt, bool print) `````` Kirill Terekhov committed Jul 18, 2016 1018 `````` { `````` Kirill Terekhov committed Jul 28, 2016 1019 `````` //Tag pnt = m->CreateTag("PROJ_PNT"+m->GetLocalProcessorRank(),DATA_REAL,NODE,NODE,2); `````` Kirill Terekhov committed Aug 12, 2016 1020 `````` std::set intersections; `````` Kirill Terekhov committed Jun 15, 2017 1021 `````` std::vector initials(segments.size() * 2); `````` Kirill Terekhov committed Jul 28, 2016 1022 `````` MarkerType initial = m->CreatePrivateMarker(); `````` Kirill Terekhov committed Feb 21, 2021 1023 `````` for (size_t k = 0; k < segments.size(); ++k) `````` Kirill Terekhov committed Jul 18, 2016 1024 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1025 1026 `````` initials[k * 2 + 0] = segments[k]->getBeg()->GetHandle(); initials[k * 2 + 1] = segments[k]->getEnd()->GetHandle(); `````` Kirill Terekhov committed Jul 28, 2016 1027 1028 `````` segments[k]->getBeg()->SetPrivateMarker(initial); segments[k]->getEnd()->SetPrivateMarker(initial); `````` Kirill Terekhov committed Jun 15, 2017 1029 1030 `````` intersections.insert(make_point(segments[k]->getBeg(), pnt)); intersections.insert(make_point(segments[k]->getEnd(), pnt)); `````` Kirill Terekhov committed Jul 18, 2016 1031 `````` } `````` Kirill Terekhov committed Feb 21, 2021 1032 `````` for (size_t i = 0; i < segments.size(); ++i) `````` Kirill Terekhov committed Jul 18, 2016 1033 `````` { `````` Kirill Terekhov committed Feb 21, 2021 1034 `````` for (size_t j = i + 1; j < segments.size(); ++j) `````` Kirill Terekhov committed Jul 18, 2016 1035 `````` { `````` Kirill Terekhov committed Jun 21, 2017 1036 `````` std::pair I = intersect_segments(m, segments[i], segments[j], intersections, pnt, print); `````` Kirill Terekhov committed Jun 15, 2017 1037 `````` if (I.first) `````` Kirill Terekhov committed Jul 18, 2016 1038 1039 `````` { ElementArray splitted_a, splitted_b; `````` Kirill Terekhov committed Jun 21, 2017 1040 1041 1042 `````` split_edges(m, Node(m, I.second.node), segments[i], segments[j], splitted_a, splitted_b, transfer, print); segments[i] = splitted_a[0]; segments[j] = splitted_b[0]; `````` Kirill Terekhov committed Jul 18, 2016 1043 1044 1045 1046 1047 `````` segments.push_back(splitted_a[1]); segments.push_back(splitted_b[1]); } } } `````` Kirill Terekhov committed Nov 03, 2016 1048 `````` //nodes.clear(); `````` Kirill Terekhov committed Jun 15, 2017 1049 `````` for (std::set::iterator it = intersections.begin(); it != intersections.end(); ++it) `````` Kirill Terekhov committed Jul 18, 2016 1050 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1051 `````` if (!m->GetPrivateMarker(it->node, initial)) `````` Kirill Terekhov committed Aug 12, 2016 1052 `````` nodes.push_back(it->node); `````` Kirill Terekhov committed Jul 18, 2016 1053 `````` } `````` Kirill Terekhov committed Jun 15, 2017 1054 1055 `````` for (int k = 0; k < (int)initials.size(); ++k) m->RemPrivateMarker(initials[k], initial); `````` Kirill Terekhov committed Jul 28, 2016 1056 `````` m->ReleasePrivateMarker(initial); `````` Kirill Terekhov committed Jul 18, 2016 1057 1058 `````` } `````` Kirill Terekhov committed Aug 12, 2016 1059 `````` void block_number_union(Element n, ElementArray & adj, Tag block, Tag write) `````` Kirill Terekhov committed Jul 08, 2016 1060 `````` { `````` Kirill Terekhov committed Aug 12, 2016 1061 `````` Storage::integer_array bn = n->IntegerArray(write); `````` Kirill Terekhov committed Jun 15, 2017 1062 `````` if (!adj.empty()) `````` Kirill Terekhov committed Jul 08, 2016 1063 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1064 1065 1066 `````` dynarray uni, tmp; //std::vector uni, tmp; for (ElementArray::size_type k = 0; k < adj.size(); ++k) `````` Kirill Terekhov committed Jul 08, 2016 1067 `````` { `````` Kirill Terekhov committed Aug 12, 2016 1068 `````` Storage::integer_array be = adj[k]->IntegerArray(block); `````` Kirill Terekhov committed Jun 15, 2017 1069 1070 `````` tmp.resize(uni.size() + be.size()); tmp.resize(std::set_union(uni.begin(), uni.end(), be.begin(), be.end(), tmp.begin()) - tmp.begin()); `````` Kirill Terekhov committed Aug 12, 2016 1071 `````` uni.swap(tmp); `````` Kirill Terekhov committed Jul 08, 2016 1072 `````` } `````` Kirill Terekhov committed Jun 15, 2017 1073 `````` bn.replace(bn.begin(), bn.end(), uni.begin(), uni.end()); `````` Kirill Terekhov committed Jul 08, 2016 1074 `````` } `````` Kirill Terekhov committed Aug 12, 2016 1075 `````` else bn.clear(); `````` Kirill Terekhov committed Jul 08, 2016 1076 `````` } `````` Kirill Terekhov committed Jun 15, 2017 1077 `````` `````` Kirill Terekhov committed Aug 12, 2016 1078 1079 1080 `````` 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 1081 `````` if (!adj.empty()) `````` Kirill Terekhov committed Jul 08, 2016 1082 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1083 1084 1085 `````` 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 1086 `````` { `````` Kirill Terekhov committed Aug 12, 2016 1087 `````` Storage::integer_array be = adj[k]->IntegerArray(block); `````` Kirill Terekhov committed Jun 15, 2017 1088 1089 `````` tmp.resize(uni.size() + be.size()); tmp.resize(std::set_union(uni.begin(), uni.end(), be.begin(), be.end(), tmp.begin()) - tmp.begin()); `````` Kirill Terekhov committed Aug 12, 2016 1090 `````` uni.swap(tmp); `````` Kirill Terekhov committed Jul 08, 2016 1091 `````` } `````` Kirill Terekhov committed Jun 15, 2017 1092 `````` bn.replace(bn.begin(), bn.end(), uni.begin(), uni.end()); `````` Kirill Terekhov committed Jul 08, 2016 1093 `````` } `````` Kirill Terekhov committed Jul 11, 2016 1094 `````` } `````` Kirill Terekhov committed Aug 12, 2016 1095 `````` `````` Kirill Terekhov committed Jun 15, 2017 1096 `````` void block_number_intersection(Element n, const ElementArray & adj, Tag block, Tag write) `````` Kirill Terekhov committed Nov 04, 2016 1097 1098 `````` { Storage::integer_array bn = n->IntegerArray(write); `````` Kirill Terekhov committed Jun 15, 2017 1099 `````` if (!adj.empty()) `````` Kirill Terekhov committed Nov 04, 2016 1100 1101 `````` { Storage::integer_array be = adj[0]->IntegerArray(block); `````` Kirill Terekhov committed Jun 15, 2017 1102 1103 1104 `````` 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 1105 1106 `````` { be = adj[k]->IntegerArray(block); `````` Kirill Terekhov committed Jun 15, 2017 1107 `````` tmp.resize(std::set_intersection(inter.begin(), inter.end(), be.begin(), be.end(), tmp.begin()) - tmp.begin()); `````` Kirill Terekhov committed Nov 04, 2016 1108 1109 `````` inter.swap(tmp); } `````` Kirill Terekhov committed Jun 15, 2017 1110 `````` bn.replace(bn.begin(), bn.end(), inter.begin(), inter.end()); `````` Kirill Terekhov committed Nov 04, 2016 1111 1112 1113 1114 1115 1116 `````` } else bn.clear(); //Storage::integer_array bn = n->IntegerArray(block); //bn.replace(bn.begin(),bn.end(),inter.begin(),inter.end()); } `````` Kirill Terekhov committed Jun 15, 2017 1117 1118 `````` struct mult_rec `````` Kirill Terekhov committed Jul 11, 2016 1119 `````` { `````` Kirill Terekhov committed Jun 15, 2017 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 `````` std::string name; double mult; int bijk[6]; }; struct editnnc_rec { int bijk[6]; double mult; }; struct fault_rec { std::string name; int dir; int bijk[6]; }; struct multflt_rec { std::string name; double mult; }; `````` Kirill Terekhov committed Jul 08, 2016 1143 `````` `````` Kirill Terekhov committed Jun 21, 2017 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