1 // 2 // Read user dependent output file 3 // 4 5 6 #include <mystdlib.h> 7 #include <myadt.hpp> 8 #include <linalg.hpp> 9 #include <csg.hpp> 10 #include <stlgeom.hpp> 11 #include <meshing.hpp> 12 13 #include "writeuser.hpp" 14 15 namespace netgen 16 { ReadFile(Mesh & mesh,const string & hfilename)17 void ReadFile (Mesh & mesh, 18 const string & hfilename) 19 { 20 PrintMessage(3, "Read User File"); 21 22 const char * filename = hfilename.c_str(); 23 24 char reco[100]; 25 int np, nbe; 26 27 28 29 // ".surf" - mesh 30 31 if ( (strlen (filename) > 5) && 32 strcmp (&filename[strlen (filename)-5], ".surf") == 0 ) 33 34 { 35 cout << "Surface file" << endl; 36 37 ifstream in (filename); 38 39 in >> reco; 40 in >> np; 41 for (int i = 1; i <= np; i++) 42 { 43 Point3d p; 44 in >> p.X() >> p.Y() >> p.Z(); 45 mesh.AddPoint (p); 46 } 47 48 mesh.ClearFaceDescriptors(); 49 mesh.AddFaceDescriptor (FaceDescriptor(1,1,0,0)); 50 51 in >> nbe; 52 // int invert = globflags.GetDefineFlag ("invertsurfacemesh"); 53 for (int i = 1; i <= nbe; i++) 54 { 55 Element2d el; 56 el.SetIndex(1); 57 58 for (int j = 1; j <= 3; j++) 59 { 60 in >> el.PNum(j); 61 // el.PNum(j)++; 62 if (el.PNum(j) < PointIndex(1) || 63 el.PNum(j) > PointIndex(np)) 64 { 65 cerr << "Point Number " << el.PNum(j) << " out of range 1..." 66 << np << endl; 67 return; 68 } 69 } 70 /* 71 if (invert) 72 swap (el.PNum(2), el.PNum(3)); 73 */ 74 75 mesh.AddSurfaceElement (el); 76 } 77 78 79 cout << "points: " << np << " faces: " << nbe << endl; 80 } 81 82 83 84 85 86 if ( (strlen (filename) > 4) && 87 strcmp (&filename[strlen (filename)-4], ".unv") == 0 ) 88 { 89 char reco[100]; 90 // int invert; 91 92 ifstream in(filename); 93 94 mesh.ClearFaceDescriptors(); 95 mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0)); 96 mesh.GetFaceDescriptor(1).SetBCProperty (1); 97 // map from unv element nr to our element number + an index if it is vol (0), bnd(1), ... 98 std::map<size_t, std::tuple<size_t, int>> element_map; 99 int dim = 3; 100 int bccounter = 0; 101 102 NgArray<Segment> tmp_segments; 103 while (in.good()) 104 { 105 in >> reco; 106 if (strcmp(reco, "-1") == 0) 107 continue; 108 109 else if (strcmp (reco, "2411") == 0) 110 { 111 cout << "nodes found" << endl; 112 113 while (1) 114 { 115 int pi, hi; 116 Point<3> p; 117 118 in >> pi; 119 if (pi == -1) 120 break; 121 122 in >> hi >> hi >> hi; 123 in >> p(0) >> p(1) >> p(2); 124 125 mesh.AddPoint (p); 126 } 127 cout << "read " << mesh.GetNP() << " points" << endl; 128 Point3d pmin, pmax; 129 mesh.GetBox (pmin, pmax); 130 if(fabs(pmin.Z() - pmax.Z()) < 1e-10 * Dist(pmin, pmax)) 131 { 132 cout << "Set Dimension to 2." << endl; 133 mesh.SetDimension(2); 134 dim = 2 ; 135 } 136 137 } 138 139 else if (strcmp (reco, "2412") == 0) 140 { 141 cout << "elements found" << endl; 142 143 while (1) 144 { 145 int label, fe_id, phys_prop, mat_prop, color, nnodes; 146 int nodes[100]; 147 int hi; 148 149 in >> label; 150 if (label == -1) break; 151 in >> fe_id >> phys_prop >> mat_prop >> color >> nnodes; 152 153 if (fe_id >= 11 && fe_id <= 32) 154 in >> hi >> hi >> hi; 155 156 157 for (int j = 0; j < nnodes; j++) 158 in >> nodes[j]; 159 160 switch (fe_id) 161 { 162 case 11: // (Rod) SEGM 163 { 164 Segment el; 165 el[0] = nodes[0]; 166 el[1] = nodes[1]; 167 el[2] = -1; 168 169 if(dim == 3){ 170 auto nr = tmp_segments.Size(); 171 tmp_segments.Append(el); 172 element_map[label] = std::make_tuple(nr+1, 2); 173 } 174 else if(dim == 2){ 175 el.si = -1; // add label to segment, will be changed later when BC's are assigned 176 auto nr = mesh.AddSegment(el); 177 element_map[label] = std::make_tuple(nr+1, 2); 178 } 179 break; 180 } 181 182 case 22: // (Tapered beam) SEGM 183 { 184 Segment el; 185 el[0] = nodes[0]; 186 el[1] = nodes[2]; 187 el[2] = nodes[1]; 188 189 if(dim == 3){ 190 auto nr = tmp_segments.Size(); 191 tmp_segments.Append(el); 192 element_map[label] = std::make_tuple(nr+1, 2); 193 } 194 else if(dim == 2){ 195 el.si = -1; // add label to segment, will be changed later when BC's are assigned 196 auto nr = mesh.AddSegment(el); 197 element_map[label] = std::make_tuple(nr+1, 2); 198 } 199 200 break; 201 } 202 case 41: // TRIG 203 { 204 Element2d el (TRIG); 205 el.SetIndex (1); 206 for (int j = 0; j < nnodes; j++) 207 el[j] = nodes[j]; 208 auto nr = mesh.AddSurfaceElement (el); 209 element_map[label] = std::make_tuple(nr+1, 1); 210 break; 211 } 212 case 42: // TRIG6 213 { 214 Element2d el(TRIG6); 215 el.SetIndex(1); 216 int jj = 0; 217 for(auto j : {0,2,4,3,5,1}) 218 el[jj++] = nodes[j]; 219 auto nr = mesh.AddSurfaceElement(el); 220 element_map[label] = std::make_tuple(nr+1, 1); 221 break; 222 } 223 case 111: // TET 224 { 225 Element el (TET); 226 el.SetIndex (1); 227 for (int j = 0; j < nnodes; j++) 228 el[j] = nodes[j]; 229 auto nr = mesh.AddVolumeElement (el); 230 element_map[label] = std::make_tuple(nr+1, 0); 231 break; 232 } 233 case 118: // TET10 234 { 235 Element el(TET10); 236 el.SetIndex(1); 237 int jj = 0; 238 for(auto j : {0,2,4,9,1,5,6,3,7,8}) 239 el[jj++] = nodes[j]; 240 auto nr = mesh.AddVolumeElement(el); 241 element_map[label] = std::make_tuple(nr+1, 0); 242 break; 243 } 244 default: 245 cout << "Do not know fe_id = " << fe_id << ", skipping it." << endl; 246 break; 247 } 248 } 249 cout << mesh.GetNE() << " elements found" << endl; 250 cout << mesh.GetNSE() << " surface elements found" << endl; 251 252 } 253 else if(strcmp (reco, "2467") == 0) 254 { 255 int matnr = 1; 256 cout << "Groups found" << endl; 257 while(in.good()) 258 { 259 int len; 260 string name; 261 in >> len; 262 if(len == -1) 263 break; 264 for(int i=0; i < 7; i++) 265 in >> len; 266 in >> name; 267 cout << len << " element are in group " << name << endl; 268 int hi, index; 269 int fdnr, ednr; 270 271 in >> hi >> index >> hi >> hi; 272 int codim = get<1>(element_map[index]); 273 // use first element to determine if boundary or volume 274 275 switch (codim) 276 { 277 case 0: 278 { 279 mesh.SetMaterial(++matnr, name); 280 mesh.VolumeElement(get<0>(element_map[index])).SetIndex(matnr); 281 break; 282 } 283 case 1: 284 { 285 if(dim == 3) 286 { 287 int bcpr = mesh.GetNFD(); 288 fdnr = mesh.AddFaceDescriptor(FaceDescriptor(bcpr, 0,0,0)); 289 mesh.GetFaceDescriptor(fdnr).SetBCProperty(bcpr+1); 290 mesh.SetBCName(bcpr, name); 291 mesh.SurfaceElement(get<0>(element_map[index])).SetIndex(fdnr); 292 bccounter++; 293 } 294 else if(dim == 2) 295 { 296 mesh.SetMaterial(matnr, name); 297 fdnr = mesh.AddFaceDescriptor(FaceDescriptor(matnr, 0,0,0)); 298 mesh.SurfaceElement(get<0>(element_map[index])).SetIndex(matnr); 299 mesh.GetFaceDescriptor(fdnr).SetBCProperty(matnr); 300 matnr++; 301 } 302 break; 303 304 } 305 case 2: 306 { 307 if(dim == 3) 308 { 309 int bcpr = mesh.GetNCD2Names()+1; 310 auto ed = EdgeDescriptor(); 311 ed.SetSurfNr(0,bcpr);//? 312 ednr = mesh.AddEdgeDescriptor(ed); 313 mesh.SetCD2Name(bcpr, name); 314 auto nr = mesh.AddSegment(tmp_segments[get<0>(element_map[index])-1]); 315 mesh[nr].edgenr = ednr+1; 316 } 317 else if(dim == 2) 318 { 319 Segment & seg = mesh.LineSegment(get<0>(element_map[index])); 320 seg.si = bccounter + 1; 321 mesh.SetBCName(bccounter, name); 322 bccounter++; 323 } 324 break; 325 326 } 327 default: 328 { 329 cout << "Codim " << codim << " not implemented yet!" << endl; 330 } 331 } 332 333 for(int i=0; i<len-1; i++) 334 { 335 in >> hi >> index >> hi >> hi; 336 switch (codim) 337 { 338 case 0: 339 mesh.VolumeElement(get<0>(element_map[index])).SetIndex(matnr); 340 break; 341 case 1: 342 if(dim == 3) mesh.SurfaceElement(get<0>(element_map[index])).SetIndex(fdnr); 343 else if (dim == 2){ 344 mesh.SurfaceElement(get<0>(element_map[index])).SetIndex(matnr-1); 345 mesh.GetFaceDescriptor(fdnr).SetBCProperty(matnr); 346 } 347 break; 348 case 2: 349 if(dim == 3) 350 { 351 auto nr = mesh.AddSegment(tmp_segments[get<0>(element_map[index])-1]); 352 mesh[nr].edgenr = ednr+1; 353 } 354 else if(dim == 2) 355 { 356 Segment & seg = mesh.LineSegment(get<0>(element_map[index])); 357 seg.si = bccounter; 358 } 359 break; 360 default: 361 break; 362 } 363 } 364 } 365 } 366 else 367 { 368 cout << "Do not know data field type " << reco << ", skipping it" << endl; 369 while(in.good()) 370 { 371 in >> reco; 372 if(strcmp(reco, "-1") == 0) 373 break; 374 } 375 } 376 } 377 378 if(dim == 2){ 379 // loop through segments to assign default BC to unmarked edges 380 int bccounter_tmp = bccounter; 381 for(int index=1; index <= mesh.GetNSeg(); index++){ 382 Segment & seg = mesh.LineSegment(index); 383 if(seg.si == -1){ 384 seg.si = bccounter + 1; 385 if(bccounter_tmp == bccounter) mesh.SetBCName(bccounter, "default"); // could be more efficient 386 bccounter_tmp++; 387 } 388 } 389 if(bccounter_tmp > bccounter) bccounter++; 390 } 391 392 393 Point3d pmin, pmax; 394 mesh.ComputeNVertices(); 395 mesh.RebuildSurfaceElementLists(); 396 mesh.GetBox (pmin, pmax); 397 mesh.UpdateTopology(); 398 if(dim == 3) bccounter++; 399 cout << "bounding-box = " << pmin << "-" << pmax << endl; 400 cout << "Created " << bccounter << " boundaries." << endl; 401 for(int i=0; i<bccounter; i++){ 402 cout << mesh.GetBCName(i) << endl; 403 } 404 } 405 406 407 408 // fepp format2d: 409 410 if ( (strlen (filename) > 7) && 411 strcmp (&filename[strlen (filename)-7], ".mesh2d") == 0 ) 412 { 413 cout << "Reading FEPP2D Mesh" << endl; 414 415 char buf[100]; 416 int np, ne, nseg, i, j; 417 418 ifstream in (filename); 419 420 in >> buf; 421 422 in >> nseg; 423 for (i = 1; i <= nseg; i++) 424 { 425 int bound, p1, p2; 426 in >> bound >> p1 >> p2; 427 // forget them 428 } 429 430 in >> ne; 431 for (i = 1; i <= ne; i++) 432 { 433 int mat, nelp; 434 in >> mat >> nelp; 435 Element2d el (nelp == 3 ? TRIG : QUAD); 436 el.SetIndex (mat); 437 for (j = 1; j <= nelp; j++) 438 in >> el.PNum(j); 439 mesh.AddSurfaceElement (el); 440 } 441 442 in >> np; 443 for (i = 1; i <= np; i++) 444 { 445 Point3d p(0,0,0); 446 in >> p.X() >> p.Y(); 447 mesh.AddPoint (p); 448 } 449 } 450 451 452 else if ( (strlen (filename) > 5) && 453 strcmp (&filename[strlen (filename)-5], ".mesh") == 0 ) 454 { 455 cout << "Reading Neutral Format" << endl; 456 457 int np, ne, nse, i, j; 458 459 ifstream in (filename); 460 461 in >> np; 462 463 if (in.good()) 464 { 465 // file starts with an integer 466 467 for (i = 1; i <= np; i++) 468 { 469 Point3d p(0,0,0); 470 in >> p.X() >> p.Y() >> p.Z(); 471 mesh.AddPoint (p); 472 } 473 474 in >> ne; 475 for (i = 1; i <= ne; i++) 476 { 477 int mat; 478 in >> mat; 479 Element el (4); 480 el.SetIndex (mat); 481 for (j = 1; j <= 4; j++) 482 in >> el.PNum(j); 483 mesh.AddVolumeElement (el); 484 } 485 486 mesh.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0)); 487 int nfd = 1; 488 489 in >> nse; 490 for (i = 1; i <= nse; i++) 491 { 492 int mat; // , nelp; 493 in >> mat; 494 Element2d el (TRIG); 495 el.SetIndex (mat); 496 while(nfd<mat) 497 { 498 ++nfd; 499 mesh.AddFaceDescriptor(FaceDescriptor(nfd,nfd,0,0)); 500 } 501 for (j = 1; j <= 3; j++) 502 in >> el.PNum(j); 503 mesh.AddSurfaceElement (el); 504 } 505 } 506 else 507 { 508 char buf[100]; 509 in.clear(); 510 do 511 { 512 in >> buf; 513 cout << "buf = " << buf << endl; 514 if (strcmp (buf, "points") == 0) 515 { 516 in >> np; 517 cout << "np = " << np << endl; 518 } 519 } 520 while (in.good()); 521 } 522 } 523 524 525 if ( (strlen (filename) > 4) && 526 strcmp (&filename[strlen (filename)-4], ".emt") == 0 ) 527 { 528 ifstream inemt (filename); 529 530 string pktfile = filename; 531 int len = strlen (filename); 532 pktfile[len-3] = 'p'; 533 pktfile[len-2] = 'k'; 534 pktfile[len-1] = 't'; 535 cout << "pktfile = " << pktfile << endl; 536 537 int np, nse, i; 538 int bcprop; 539 ifstream inpkt (pktfile.c_str()); 540 inpkt >> np; 541 NgArray<double> values(np); 542 for (i = 1; i <= np; i++) 543 { 544 Point3d p(0,0,0); 545 inpkt >> p.X() >> p.Y() >> p.Z() 546 >> bcprop >> values.Elem(i); 547 mesh.AddPoint (p); 548 } 549 550 mesh.ClearFaceDescriptors(); 551 mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0)); 552 mesh.GetFaceDescriptor(1).SetBCProperty (1); 553 mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0)); 554 mesh.GetFaceDescriptor(2).SetBCProperty (2); 555 mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0)); 556 mesh.GetFaceDescriptor(3).SetBCProperty (3); 557 mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0)); 558 mesh.GetFaceDescriptor(4).SetBCProperty (4); 559 mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0)); 560 mesh.GetFaceDescriptor(5).SetBCProperty (5); 561 562 int p1, p2, p3; 563 double value; 564 inemt >> nse; 565 for (i = 1; i <= nse; i++) 566 { 567 inemt >> p1 >> p2 >> p3 >> bcprop >> value; 568 569 if (bcprop < 1 || bcprop > 4) 570 cerr << "bcprop out of range, bcprop = " << bcprop << endl; 571 p1++; 572 p2++; 573 p3++; 574 if (p1 < 1 || p1 > np || p2 < 1 || p2 > np || p3 < 1 || p3 > np) 575 { 576 cout << "p1 = " << p1 << " p2 = " << p2 << " p3 = " << p3 << endl; 577 } 578 579 if (i > 110354) Swap (p2, p3); 580 if (mesh.Point(p1)(0) < 0.25) 581 Swap (p2,p3); 582 583 Element2d el(TRIG); 584 585 if (bcprop == 1) 586 { 587 if (values.Get(p1) < -69999) 588 el.SetIndex(1); 589 else 590 el.SetIndex(2); 591 } 592 else 593 el.SetIndex(3); 594 595 596 el.PNum(1) = p1; 597 el.PNum(2) = p2; 598 el.PNum(3) = p3; 599 mesh.AddSurfaceElement (el); 600 } 601 602 603 ifstream incyl ("ngusers/guenter/cylinder.surf"); 604 int npcyl, nsecyl; 605 incyl >> npcyl; 606 cout << "npcyl = " << npcyl << endl; 607 for (i = 1; i <= npcyl; i++) 608 { 609 Point3d p(0,0,0); 610 incyl >> p.X() >> p.Y() >> p.Z(); 611 mesh.AddPoint (p); 612 } 613 incyl >> nsecyl; 614 cout << "nsecyl = " << nsecyl << endl; 615 for (i = 1; i <= nsecyl; i++) 616 { 617 incyl >> p1 >> p2 >> p3; 618 p1 += np; 619 p2 += np; 620 p3 += np; 621 Element2d el(TRIG); 622 el.SetIndex(5); 623 el.PNum(1) = p1; 624 el.PNum(2) = p2; 625 el.PNum(3) = p3; 626 mesh.AddSurfaceElement (el); 627 } 628 } 629 630 631 // .tet mesh 632 if ( (strlen (filename) > 4) && 633 strcmp (&filename[strlen (filename)-4], ".tet") == 0 ) 634 { 635 ReadTETFormat (mesh, filename); 636 } 637 638 639 // .fnf mesh (FNF - PE neutral format) 640 if ( (strlen (filename) > 4) && 641 strcmp (&filename[strlen (filename)-4], ".fnf") == 0 ) 642 { 643 ReadFNFFormat (mesh, filename); 644 } 645 646 // .cgns file - CFD General Notation System 647 if ( (strlen (filename) > 5) && 648 strcmp (&filename[strlen (filename)-5], ".cgns") == 0 ) 649 { 650 ReadCGNSMesh (mesh, filename); 651 } 652 653 if ( ( (strlen (filename) > 4) && strcmp (&filename[strlen (filename)-4], ".stl") == 0 ) || 654 ( (strlen (filename) > 5) && strcmp (&filename[strlen (filename)-5], ".stlb") == 0 ) ) 655 { 656 ifstream ist{string{filename}}; 657 auto geom = shared_ptr<STLGeometry>(STLGeometry::Load(ist)); 658 659 mesh.SetDimension (3); 660 661 auto & points = geom->GetPoints(); 662 663 for (auto & p : points) 664 mesh.AddPoint(MeshPoint(p)); 665 666 mesh.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 1)); 667 668 for (auto ti : IntRange(geom->GetNT())) 669 { 670 Element2d el(TRIG); 671 for (auto i : IntRange(3)) 672 el[i] = int((*geom)[STLTrigId(ti+IndexBASE<netgen::STLTrigId>())][i]); 673 674 el.SetIndex(1); 675 676 mesh.AddSurfaceElement(el); 677 } 678 } 679 } 680 681 } 682 683