1 // Copyright(C) 1999-2020 National Technology & Engineering Solutions 2 // of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with 3 // NTESS, the U.S. Government retains certain rights in this software. 4 // 5 // See packages/seacas/LICENSE for details 6 7 #include <Ioss_EntityType.h> // for EntityType, etc 8 #include <Ioss_Hex8.h> 9 #include <Ioss_Utils.h> 10 #include <algorithm> 11 #include <cassert> // for assert 12 #include <cmath> // for atan2, cos, sin 13 #include <cstdlib> // for nullptr, exit, etc 14 #include <fmt/ostream.h> 15 #include <gen_struc/Iogs_GeneratedMesh.h> 16 #include <numeric> 17 #include <string> 18 #include <sys/types.h> // for ssize_t 19 #include <tokenize.h> // for tokenize 20 #include <vector> // for vector 21 22 namespace Iogs { GeneratedMesh(int64_t,int64_t,int64_t,int proc_count,int my_proc)23 GeneratedMesh::GeneratedMesh(int64_t /*num_x */, int64_t /* num_y */, int64_t /* num_z */, 24 int proc_count, int my_proc) 25 : processorCount(proc_count), myProcessor(my_proc) 26 { 27 initialize(); 28 } 29 GeneratedMesh(const std::string & parameters,int proc_count,int my_proc)30 GeneratedMesh::GeneratedMesh(const std::string ¶meters, int proc_count, int my_proc) 31 : processorCount(proc_count), myProcessor(my_proc) 32 { 33 // Possible that the 'parameters' has the working directory path 34 // prepended to the parameter list. Strip off everything in front 35 // of the last '/' (if any)... 36 auto params = Ioss::tokenize(parameters, "/"); 37 38 auto groups = Ioss::tokenize(params.back(), "|+"); 39 40 // First 'group' is the interval specification -- IxJxK 41 auto tokens = Ioss::tokenize(groups[0], "x"); 42 assert(tokens.size() == 3); 43 numX = std::stoll(tokens[0]); 44 numY = std::stoll(tokens[1]); 45 numZ = std::stoll(tokens[2]); 46 47 initialize(); 48 parse_options(groups); 49 } 50 GeneratedMesh()51 GeneratedMesh::GeneratedMesh() { initialize(); } 52 53 GeneratedMesh::~GeneratedMesh() = default; 54 initialize()55 void GeneratedMesh::initialize() 56 { 57 if (processorCount > numZ) { 58 std::ostringstream errmsg; 59 fmt::print(errmsg, 60 "ERROR: ({})\n" 61 " The number of mesh intervals in the Z direction ({})\n" 62 " must be at least as large as the number of processors ({}).\n" 63 " The current parameters do not meet that requirement. Execution will " 64 "terminate.\n", 65 __func__, numZ, processorCount); 66 IOSS_ERROR(errmsg); 67 } 68 69 if (processorCount > 1) { 70 myNumZ = numZ / processorCount; 71 if (myProcessor < (numZ % processorCount)) { 72 myNumZ++; 73 } 74 75 // Determine myStartZ for this processor... 76 size_t extra = numZ % processorCount; 77 if (extra > myProcessor) { 78 extra = myProcessor; 79 } 80 size_t per_proc = numZ / processorCount; 81 myStartZ = myProcessor * per_proc + extra; 82 } 83 else { 84 myNumZ = numZ; 85 } 86 87 for (int i = 0; i < 3; i++) { 88 for (int j = 0; j < 3; j++) { 89 rotmat[i][j] = 0.0; 90 } 91 rotmat[i][i] = 1.0; 92 } 93 94 variableCount[Ioss::COMMSET] = 0; 95 variableCount[Ioss::EDGEBLOCK] = 0; 96 variableCount[Ioss::EDGESET] = 0; 97 variableCount[Ioss::ELEMENTBLOCK] = 0; 98 variableCount[Ioss::ELEMENTSET] = 0; 99 variableCount[Ioss::FACEBLOCK] = 0; 100 variableCount[Ioss::FACESET] = 0; 101 variableCount[Ioss::INVALID_TYPE] = 0; 102 variableCount[Ioss::NODEBLOCK] = 0; 103 variableCount[Ioss::REGION] = 0; 104 variableCount[Ioss::SIDEBLOCK] = 0; 105 variableCount[Ioss::SIDESET] = 0; 106 variableCount[Ioss::SUPERELEMENT] = 0; 107 } 108 add_sideset(ShellLocation loc)109 int64_t GeneratedMesh::add_sideset(ShellLocation loc) 110 { 111 sidesets.push_back(loc); 112 return sidesets.size(); 113 } 114 set_bbox(double xmin,double ymin,double zmin,double xmax,double ymax,double zmax)115 void GeneratedMesh::set_bbox(double xmin, double ymin, double zmin, double xmax, double ymax, 116 double zmax) 117 { 118 // NOTE: All calculations are based on the currently 119 // active interval settings. If scale or offset or zdecomp 120 // specified later in the option list, you may not get the 121 // desired bounding box. 122 if (numX == 0 || numY == 0 || numZ == 0) { 123 std::ostringstream errmsg; 124 fmt::print(errmsg, 125 "ERROR: ({})\n" 126 " All interval counts must be greater than 0.\n" 127 " numX = {}, numY = {}, numZ = {}\n", 128 __func__, numX, numY, numZ); 129 IOSS_ERROR(errmsg); 130 } 131 132 double x_range = xmax - xmin; 133 double y_range = ymax - ymin; 134 double z_range = zmax - zmin; 135 136 sclX = x_range / static_cast<double>(numX); 137 sclY = y_range / static_cast<double>(numY); 138 sclZ = z_range / static_cast<double>(numZ); 139 140 offX = xmin; 141 offY = ymin; 142 offZ = zmin; 143 } 144 set_scale(double scl_x,double scl_y,double scl_z)145 void GeneratedMesh::set_scale(double scl_x, double scl_y, double scl_z) 146 { 147 sclX = scl_x; 148 sclY = scl_y; 149 sclZ = scl_z; 150 } 151 set_offset(double off_x,double off_y,double off_z)152 void GeneratedMesh::set_offset(double off_x, double off_y, double off_z) 153 { 154 offX = off_x; 155 offY = off_y; 156 offZ = off_z; 157 } 158 parse_options(const std::vector<std::string> & groups)159 void GeneratedMesh::parse_options(const std::vector<std::string> &groups) 160 { 161 for (size_t i = 1; i < groups.size(); i++) { 162 auto option = Ioss::tokenize(groups[i], ":"); 163 // option[0] is the type of the option and option[1] is the argument to the option. 164 165 if (option[0] == "sideset") { 166 // Option of the form "sideset:xXyYzZ" 167 // The argument specifies whether there is a sideset 168 // at the location. 'x' is minX, 'X' is maxX, etc. 169 for (auto opt : option[1]) { 170 switch (opt) { 171 case 'x': add_sideset(MX); break; 172 case 'X': add_sideset(PX); break; 173 case 'y': add_sideset(MY); break; 174 case 'Y': add_sideset(PY); break; 175 case 'z': add_sideset(MZ); break; 176 case 'Z': add_sideset(PZ); break; 177 default: 178 std::ostringstream errmsg; 179 fmt::print(errmsg, "ERROR: Unrecognized sideset location option '{}'.", opt); 180 IOSS_ERROR(errmsg); 181 } 182 } 183 } 184 else if (option[0] == "scale") { 185 // Option of the form "scale:xs,ys,zs 186 auto tokens = Ioss::tokenize(option[1], ","); 187 assert(tokens.size() == 3); 188 sclX = std::stod(tokens[0]); 189 sclY = std::stod(tokens[1]); 190 sclZ = std::stod(tokens[2]); 191 } 192 193 else if (option[0] == "offset") { 194 // Option of the form "offset:xo,yo,zo 195 auto tokens = Ioss::tokenize(option[1], ","); 196 assert(tokens.size() == 3); 197 offX = std::stod(tokens[0]); 198 offY = std::stod(tokens[1]); 199 offZ = std::stod(tokens[2]); 200 } 201 202 else if (option[0] == "zdecomp") { 203 // Option of the form "zdecomp:1,1,2,2,1,2,... 204 // Specifies the number of intervals in the z direction 205 // for each processor. The number of tokens must match 206 // the number of processors. Note that the new numZ will 207 // be the sum of the intervals specified in this command. 208 auto tokens = Ioss::tokenize(option[1], ","); 209 assert(tokens.size() == processorCount); 210 Ioss::Int64Vector Zs; 211 numZ = 0; 212 for (size_t j = 0; j < processorCount; j++) { 213 Zs.push_back(std::stoll(tokens[j])); 214 numZ += Zs[j]; 215 } 216 myNumZ = Zs[myProcessor]; 217 myStartZ = 0; 218 for (size_t j = 0; j < myProcessor; j++) { 219 myStartZ += Zs[j]; 220 } 221 } 222 223 else if (option[0] == "bbox") { 224 // Bounding-Box Option of the form "bbox:xmin,ymin,zmin,xmax,ymax,zmaxo 225 auto tokens = Ioss::tokenize(option[1], ","); 226 assert(tokens.size() == 6); 227 double xmin = std::stod(tokens[0]); 228 double ymin = std::stod(tokens[1]); 229 double zmin = std::stod(tokens[2]); 230 double xmax = std::stod(tokens[3]); 231 double ymax = std::stod(tokens[4]); 232 double zmax = std::stod(tokens[5]); 233 234 set_bbox(xmin, ymin, zmin, xmax, ymax, zmax); 235 } 236 237 else if (option[0] == "rotate") { 238 // Rotate Option of the form "rotate:axis,angle,axis,angle,... 239 auto tokens = Ioss::tokenize(option[1], ","); 240 assert(tokens.size() % 2 == 0); 241 for (size_t ir = 0; ir < tokens.size();) { 242 std::string axis = tokens[ir++]; 243 double angle_degree = std::stod(tokens[ir++]); 244 set_rotation(axis, angle_degree); 245 } 246 } 247 248 else if (option[0] == "times") { 249 timestepCount = std::stoll(option[1]); 250 } 251 252 else if (option[0] == "variables") { 253 // Variables Option of the form "variables:global,10,element,100,..." 254 auto tokens = Ioss::tokenize(option[1], ","); 255 assert(tokens.size() % 2 == 0); 256 for (size_t ir = 0; ir < tokens.size();) { 257 std::string type = tokens[ir++]; 258 int count = std::stoi(tokens[ir++]); 259 set_variable_count(type, count); 260 } 261 if (timestepCount == 0) { 262 timestepCount = 1; 263 } 264 } 265 266 else if (option[0] == "help") { 267 fmt::print(Ioss::OUTPUT(), 268 "\nValid Options for GeneratedMesh parameter string:\n" 269 "\tIxJxK -- specifies intervals; must be first option. Ex: 4x10x12\n" 270 "\toffset:xoff, yoff, zoff\n" 271 "\tscale: xscl, yscl, zscl\n" 272 "\tzdecomp:n1,n2,n3,...,n#proc\n" 273 "\tbbox: xmin, ymin, zmin, xmax, ymax, zmax\n" 274 "\trotate: axis,angle,axis,angle,...\n" 275 "\tsideset:xXyYzZ (specifies which plane to apply sideset)\n" 276 "\tvariables:type,count,... " 277 "type=global|element|node|nodal|sideset|surface\n" 278 "\ttimes:count (number of timesteps to generate)\n" 279 "\tshow -- show mesh parameters\n" 280 "\thelp -- show this list\n\n"); 281 } 282 283 else if (option[0] == "show") { 284 show_parameters(); 285 } 286 else { 287 fmt::print(Ioss::WARNING(), "Unrecognized option '{}'. It will be ignored.\n", option[0]); 288 } 289 } 290 } 291 show_parameters()292 void GeneratedMesh::show_parameters() const 293 { 294 if (myProcessor == 0) { 295 fmt::print(Ioss::OUTPUT(), 296 "\nMesh Parameters:\n" 297 "\tIntervals: {} by {} by {}\n" 298 "\tX = {} * (0..{}) + {}\tRange: {} <= X <= {}\n" 299 "\tY = {} * (0..{}) + {}\tRange: {} <= Y <= {}\n" 300 "\tZ = {} * (0..{}) + {}\tRange: {} <= Z <= {}\n\n" 301 "\tNode Count (total) = {:12L}\n" 302 "\tCell Count (total) = {:12L}\n" 303 "\tBlock Count = {:12L}\n" 304 "\tSideSet Count = {:12L}\n" 305 "\tTimestep Count = {:12L}\n\n", 306 numX, numY, numZ, sclX, numX, offX, offX, offX + numX * sclX, sclY, numY, offY, 307 offY, offY + numY * sclY, sclZ, numZ, offZ, offZ, offZ + numZ * sclZ, node_count(), 308 element_count(), structured_block_count(), sideset_count(), timestep_count()); 309 310 if (doRotation) { 311 fmt::print(Ioss::OUTPUT(), "\tRotation Matrix: \n\t"); 312 for (auto &elem : rotmat) { 313 for (double jj : elem) { 314 fmt::print(Ioss::OUTPUT(), "{:14.e}\t", jj); 315 } 316 fmt::print(Ioss::OUTPUT(), "\n\t"); 317 } 318 fmt::print(Ioss::OUTPUT(), "\n"); 319 } 320 } 321 } 322 node_count()323 int64_t GeneratedMesh::node_count() const { return (numX + 1) * (numY + 1) * (numZ + 1); } 324 node_count_proc()325 int64_t GeneratedMesh::node_count_proc() const { return (numX + 1) * (numY + 1) * (myNumZ + 1); } 326 structured_block_count()327 int64_t GeneratedMesh::structured_block_count() const { return 1; } 328 sideset_count()329 int64_t GeneratedMesh::sideset_count() const { return sidesets.size(); } 330 element_count()331 int64_t GeneratedMesh::element_count() const 332 { 333 int64_t count = element_count(1); 334 return count; 335 } 336 element_count_proc()337 int64_t GeneratedMesh::element_count_proc() const 338 { 339 int64_t count = 0; 340 for (int64_t i = 0; i < structured_block_count(); i++) { 341 count += element_count_proc(i + 1); 342 } 343 return count; 344 } 345 element_count(int64_t block_number)346 int64_t GeneratedMesh::element_count(int64_t block_number) const 347 { 348 assert(block_number <= structured_block_count()); 349 return numX * numY * numZ; 350 } 351 element_count_proc(int64_t block_number)352 int64_t GeneratedMesh::element_count_proc(int64_t block_number) const 353 { 354 assert(block_number <= structured_block_count()); 355 return numX * numY * myNumZ; 356 } 357 sideset_side_count(int64_t id)358 int64_t GeneratedMesh::sideset_side_count(int64_t id) const 359 { 360 // id is position in sideset list + 1 361 assert(id > 0 && (size_t)id <= sidesets.size()); 362 ShellLocation loc = sidesets[id - 1]; 363 switch (loc) { 364 case MX: 365 case PX: return numY * numZ; 366 case MY: 367 case PY: return numX * numZ; 368 case MZ: 369 case PZ: return numX * numY; 370 } 371 return 0; 372 } 373 sideset_side_count_proc(int64_t id)374 int64_t GeneratedMesh::sideset_side_count_proc(int64_t id) const 375 { 376 // id is position in sideset list + 1 377 assert(id > 0 && (size_t)id <= sidesets.size()); 378 ShellLocation loc = sidesets[id - 1]; 379 switch (loc) { 380 case MX: 381 case PX: return numY * myNumZ; 382 case MY: 383 case PY: return numX * myNumZ; 384 case MZ: 385 if (myProcessor == 0) { 386 return numX * numY; 387 } 388 else { 389 return 0; 390 } 391 case PZ: 392 if (myProcessor == processorCount - 1) { 393 return numX * numY; 394 } 395 else { 396 return 0; 397 } 398 } 399 return 0; 400 } 401 topology_type(int64_t block_number)402 std::pair<std::string, int> GeneratedMesh::topology_type(int64_t block_number) const 403 { 404 assert(block_number <= structured_block_count() && block_number > 0); 405 return std::make_pair(std::string(Ioss::Hex8::name), 8); 406 } 407 node_map(Ioss::Int64Vector & map)408 void GeneratedMesh::node_map(Ioss::Int64Vector &map) const 409 { 410 map.resize(node_count_proc()); 411 int64_t offset = myStartZ * (numX + 1) * (numY + 1); 412 std::iota(map.begin(), map.end(), offset + 1); 413 } 414 node_map(Ioss::IntVector & map)415 void GeneratedMesh::node_map(Ioss::IntVector &map) const 416 { 417 map.resize(node_count_proc()); 418 int offset = myStartZ * (numX + 1) * (numY + 1); 419 std::iota(map.begin(), map.end(), offset + 1); 420 } 421 communication_node_count_proc()422 int64_t GeneratedMesh::communication_node_count_proc() const 423 { 424 int64_t count = (numX + 1) * (numY + 1); 425 if (myProcessor != 0 && myProcessor != processorCount - 1) { 426 count *= 2; 427 } 428 429 return count; 430 } 431 owning_processor(int * owner,int64_t num_node)432 void GeneratedMesh::owning_processor(int *owner, int64_t num_node) 433 { 434 for (int64_t i = 0; i < num_node; i++) { 435 owner[i] = myProcessor; 436 } 437 438 if (myProcessor != 0) { 439 int64_t count = (numX + 1) * (numY + 1); 440 for (int64_t i = 0; i < count; i++) { 441 owner[i] = myProcessor - 1; 442 } 443 } 444 } 445 build_node_map(Ioss::Int64Vector & map,std::vector<int> & proc,int64_t slab,size_t slabOffset,size_t adjacentProc,size_t index)446 void GeneratedMesh::build_node_map(Ioss::Int64Vector &map, std::vector<int> &proc, int64_t slab, 447 size_t slabOffset, size_t adjacentProc, size_t index) 448 { 449 int64_t offset = (myStartZ + slabOffset) * (numX + 1) * (numY + 1); 450 for (int64_t i = 0; i < slab; i++) { 451 map[index] = offset + i + 1; 452 proc[index++] = static_cast<int>(adjacentProc); 453 } 454 } 455 node_communication_map(Ioss::Int64Vector & map,std::vector<int> & proc)456 void GeneratedMesh::node_communication_map(Ioss::Int64Vector &map, std::vector<int> &proc) 457 { 458 bool isFirstProc = myProcessor == 0; 459 bool isLastProc = myProcessor == processorCount - 1; 460 461 int64_t count = (numX + 1) * (numY + 1); 462 int64_t slab = count; 463 if (!isFirstProc && !isLastProc) { 464 count *= 2; 465 } 466 map.resize(count); 467 proc.resize(count); 468 469 size_t offset = 0; 470 if (!isFirstProc) { 471 build_node_map(map, proc, slab, 0, myProcessor - 1, offset); 472 offset += slab; 473 } 474 if (!isLastProc) { 475 build_node_map(map, proc, slab, myNumZ, myProcessor + 1, offset); 476 } 477 } 478 element_map(int64_t block_number,Ioss::Int64Vector & map)479 void GeneratedMesh::element_map(int64_t block_number, Ioss::Int64Vector &map) const 480 { 481 raw_element_map(block_number, map); 482 } 483 element_map(int64_t block_number,Ioss::IntVector & map)484 void GeneratedMesh::element_map(int64_t block_number, Ioss::IntVector &map) const 485 { 486 raw_element_map(block_number, map); 487 } 488 489 template <typename INT> raw_element_map(int64_t block_number,std::vector<INT> & map)490 void GeneratedMesh::raw_element_map(int64_t block_number, std::vector<INT> &map) const 491 { 492 assert(block_number <= structured_block_count() && block_number > 0); 493 494 INT count = element_count_proc(block_number); 495 map.reserve(count); 496 497 if (block_number == 1) { 498 // Hex/Tet block... 499 count = element_count_proc(1); 500 INT offset = myStartZ * numX * numY; 501 for (INT i = 0; i < count; i++) { 502 map.push_back(offset + i + 1); 503 } 504 } 505 } 506 element_map(Ioss::Int64Vector & map)507 void GeneratedMesh::element_map(Ioss::Int64Vector &map) const { raw_element_map(map); } 508 element_map(Ioss::IntVector & map)509 void GeneratedMesh::element_map(Ioss::IntVector &map) const { raw_element_map(map); } 510 raw_element_map(std::vector<INT> & map)511 template <typename INT> void GeneratedMesh::raw_element_map(std::vector<INT> &map) const 512 { 513 INT count = element_count_proc(); 514 map.reserve(count); 515 516 // Hex block... 517 count = element_count_proc(1); 518 INT offset = myStartZ * numX * numY; 519 for (INT i = 0; i < count; i++) { 520 map.push_back(offset + i + 1); 521 } 522 } 523 element_surface_map(ShellLocation loc,Ioss::Int64Vector & map)524 void GeneratedMesh::element_surface_map(ShellLocation loc, Ioss::Int64Vector &map) const 525 { 526 int64_t count = 0; 527 map.resize(2 * count); 528 int64_t index = 0; 529 int64_t offset = 0; 530 531 // For hex elements 532 switch (loc) { 533 case MX: 534 offset = myStartZ * numX * numY + 1; // 1-based elem id 535 for (size_t k = 0; k < myNumZ; ++k) { 536 for (size_t j = 0; j < numY; ++j) { 537 map[index++] = offset; 538 map[index++] = 3; // 0-based local face id 539 offset += numX; 540 } 541 } 542 break; 543 544 case PX: 545 offset = myStartZ * numX * numY + numX; 546 for (size_t k = 0; k < myNumZ; ++k) { 547 for (size_t j = 0; j < numY; ++j) { 548 map[index++] = offset; // 1-based elem id 549 map[index++] = 1; // 0-based local face id 550 offset += numX; 551 } 552 } 553 break; 554 555 case MY: 556 offset = myStartZ * numX * numY + 1; 557 for (size_t k = 0; k < myNumZ; ++k) { 558 for (size_t i = 0; i < numX; ++i) { 559 map[index++] = offset++; 560 map[index++] = 0; // 0-based local face id 561 } 562 offset += numX * (numY - 1); 563 } 564 break; 565 566 case PY: 567 offset = myStartZ * numX * numY + numX * (numY - 1) + 1; 568 for (size_t k = 0; k < myNumZ; ++k) { 569 for (size_t i = 0; i < numX; ++i) { 570 map[index++] = offset++; 571 map[index++] = 2; // 0-based local face id 572 } 573 offset += numX * (numY - 1); 574 } 575 break; 576 577 case MZ: 578 if (myProcessor == 0) { 579 offset = 1; 580 for (size_t i = 0; i < numY; i++) { 581 for (size_t j = 0; j < numX; j++) { 582 map[index++] = offset++; 583 map[index++] = 4; 584 } 585 } 586 } 587 break; 588 589 case PZ: 590 if (myProcessor == processorCount - 1) { 591 offset = (numZ - 1) * numX * numY + 1; 592 for (size_t i = 0, k = 0; i < numY; i++) { 593 for (size_t j = 0; j < numX; j++, k++) { 594 map[index++] = offset++; 595 map[index++] = 5; 596 } 597 } 598 } 599 break; 600 } 601 } 602 coordinates(std::vector<double> & coord)603 void GeneratedMesh::coordinates(std::vector<double> &coord) const 604 { 605 /* create global coordinates */ 606 int64_t count = node_count_proc(); 607 coord.resize(count * 3); 608 coordinates(&coord[0]); 609 } 610 coordinates(double * coord)611 void GeneratedMesh::coordinates(double *coord) const 612 { 613 /* create global coordinates */ 614 int64_t count = node_count_proc(); 615 616 int64_t k = 0; 617 for (size_t m = myStartZ; m < myStartZ + myNumZ + 1; m++) { 618 for (size_t i = 0; i < numY + 1; i++) { 619 for (size_t j = 0; j < numX + 1; j++) { 620 coord[k++] = sclX * static_cast<double>(j) + offX; 621 coord[k++] = sclY * static_cast<double>(i) + offY; 622 coord[k++] = sclZ * static_cast<double>(m) + offZ; 623 } 624 } 625 } 626 627 if (doRotation) { 628 for (int64_t i = 0; i < count * 3; i += 3) { 629 double xn = coord[i + 0]; 630 double yn = coord[i + 1]; 631 double zn = coord[i + 2]; 632 coord[i + 0] = xn * rotmat[0][0] + yn * rotmat[1][0] + zn * rotmat[2][0]; 633 coord[i + 1] = xn * rotmat[0][1] + yn * rotmat[1][1] + zn * rotmat[2][1]; 634 coord[i + 2] = xn * rotmat[0][2] + yn * rotmat[1][2] + zn * rotmat[2][2]; 635 } 636 } 637 } 638 coordinates(std::vector<double> & x,std::vector<double> & y,std::vector<double> & z)639 void GeneratedMesh::coordinates(std::vector<double> &x, std::vector<double> &y, 640 std::vector<double> &z) const 641 { 642 /* create global coordinates */ 643 int64_t count = node_count_proc(); 644 x.reserve(count); 645 y.reserve(count); 646 z.reserve(count); 647 648 for (size_t m = myStartZ; m < myStartZ + myNumZ + 1; m++) { 649 for (size_t i = 0; i < numY + 1; i++) { 650 for (size_t j = 0; j < numX + 1; j++) { 651 x.push_back(sclX * static_cast<double>(j) + offX); 652 y.push_back(sclY * static_cast<double>(i) + offY); 653 z.push_back(sclZ * static_cast<double>(m) + offZ); 654 } 655 } 656 } 657 if (doRotation) { 658 for (int64_t i = 0; i < count; i++) { 659 double xn = x[i]; 660 double yn = y[i]; 661 double zn = z[i]; 662 x.push_back(xn * rotmat[0][0] + yn * rotmat[1][0] + zn * rotmat[2][0]); 663 y.push_back(xn * rotmat[0][1] + yn * rotmat[1][1] + zn * rotmat[2][1]); 664 z.push_back(xn * rotmat[0][2] + yn * rotmat[1][2] + zn * rotmat[2][2]); 665 } 666 } 667 } 668 coordinates(int component,std::vector<double> & xyz)669 void GeneratedMesh::coordinates(int component, std::vector<double> &xyz) const 670 { 671 assert(!doRotation); 672 /* create global coordinates */ 673 size_t count = node_count_proc(); 674 xyz.reserve(count); 675 676 double offset = 0; 677 double scale = 1; 678 if (component == 1) { 679 offset = offX; 680 scale = sclX; 681 for (size_t m = myStartZ; m < myStartZ + myNumZ + 1; m++) { 682 for (size_t i = 0; i < numY + 1; i++) { 683 for (size_t j = 0; j < numX + 1; j++) { 684 xyz.push_back(scale * static_cast<double>(j) + offset); 685 } 686 } 687 } 688 } 689 else if (component == 2) { 690 offset = offY; 691 scale = sclY; 692 for (size_t m = myStartZ; m < myStartZ + myNumZ + 1; m++) { 693 for (size_t i = 0; i < numY + 1; i++) { 694 for (size_t j = 0; j < numX + 1; j++) { 695 xyz.push_back(scale * static_cast<double>(i) + offset); 696 } 697 } 698 } 699 } 700 else if (component == 3) { 701 offset = offZ; 702 scale = sclZ; 703 for (size_t m = myStartZ; m < myStartZ + myNumZ + 1; m++) { 704 for (size_t i = 0; i < numY + 1; i++) { 705 for (size_t j = 0; j < numX + 1; j++) { 706 xyz.push_back(scale * static_cast<double>(m) + offset); 707 } 708 } 709 } 710 } 711 } 712 coordinates(int component,int,double * xyz)713 void GeneratedMesh::coordinates(int component, int /* zone */, double *xyz) const 714 { 715 assert(!doRotation); 716 /* create global coordinates */ 717 if (component == 0) { 718 size_t jjj = 0; 719 for (size_t m = 0; m < numZ + 1; m++) { 720 for (size_t i = 0; i < numY + 1; i++) { 721 for (size_t j = 0; j < numX + 1; j++) { 722 xyz[jjj++] = (sclX * static_cast<double>(j) + offX); 723 xyz[jjj++] = (sclY * static_cast<double>(i) + offY); 724 xyz[jjj++] = (sclZ * static_cast<double>(m) + offZ); 725 } 726 } 727 } 728 } 729 else if (component == 1) { 730 size_t jjj = 0; 731 for (size_t m = 0; m < numZ + 1; m++) { 732 for (size_t i = 0; i < numY + 1; i++) { 733 for (size_t j = 0; j < numX + 1; j++) { 734 xyz[jjj++] = (sclX * static_cast<double>(j) + offX); 735 } 736 } 737 } 738 } 739 else if (component == 2) { 740 size_t jjj = 0; 741 for (size_t m = 0; m < numZ + 1; m++) { 742 for (size_t i = 0; i < numY + 1; i++) { 743 for (size_t j = 0; j < numX + 1; j++) { 744 xyz[jjj++] = (sclY * static_cast<double>(i) + offY); 745 } 746 } 747 } 748 } 749 else if (component == 3) { 750 size_t jjj = 0; 751 for (size_t m = 0; m < numZ + 1; m++) { 752 for (size_t i = 0; i < numY + 1; i++) { 753 for (size_t j = 0; j < numX + 1; j++) { 754 xyz[jjj++] = (sclZ * static_cast<double>(m) + offZ); 755 } 756 } 757 } 758 } 759 } 760 connectivity(int64_t block_number,Ioss::Int64Vector & connect)761 void GeneratedMesh::connectivity(int64_t block_number, Ioss::Int64Vector &connect) const 762 { 763 if (block_number == 1) { // HEX Element Block 764 connect.resize(element_count_proc(block_number) * 8); 765 } 766 raw_connectivity(block_number, &connect[0]); 767 } 768 connectivity(int64_t block_number,Ioss::IntVector & connect)769 void GeneratedMesh::connectivity(int64_t block_number, Ioss::IntVector &connect) const 770 { 771 if (block_number == 1) { // HEX Element Block 772 connect.resize(element_count_proc(block_number) * 8); 773 } 774 raw_connectivity(block_number, &connect[0]); 775 } 776 connectivity(int64_t block_number,int64_t * connect)777 void GeneratedMesh::connectivity(int64_t block_number, int64_t *connect) const 778 { 779 raw_connectivity(block_number, connect); 780 } 781 connectivity(int64_t block_number,int * connect)782 void GeneratedMesh::connectivity(int64_t block_number, int *connect) const 783 { 784 raw_connectivity(block_number, connect); 785 } 786 787 template <typename INT> raw_connectivity(int64_t block_number,INT * connect)788 void GeneratedMesh::raw_connectivity(int64_t block_number, INT *connect) const 789 { 790 assert(block_number <= structured_block_count()); 791 792 INT xp1yp1 = (numX + 1) * (numY + 1); 793 794 /* build connectivity array (node list) for mesh */ 795 if (block_number == 1) { // main block elements 796 797 // Hex elements 798 size_t cnt = 0; 799 for (size_t m = myStartZ; m < myNumZ + myStartZ; m++) { 800 for (size_t i = 0, k = 0; i < numY; i++) { 801 for (size_t j = 0; j < numX; j++, k++) { 802 size_t base = (m * xp1yp1) + k + i + 1; 803 804 connect[cnt++] = base; 805 connect[cnt++] = base + 1; 806 connect[cnt++] = base + numX + 2; 807 connect[cnt++] = base + numX + 1; 808 809 connect[cnt++] = xp1yp1 + base; 810 connect[cnt++] = xp1yp1 + base + 1; 811 connect[cnt++] = xp1yp1 + base + numX + 2; 812 connect[cnt++] = xp1yp1 + base + numX + 1; 813 } 814 } 815 } 816 } 817 return; 818 } 819 sideset_elem_sides(int64_t id,Ioss::Int64Vector & elem_sides)820 void GeneratedMesh::sideset_elem_sides(int64_t id, Ioss::Int64Vector &elem_sides) const 821 { 822 // id is position in sideset list + 1 823 assert(id > 0 && (size_t)id <= sidesets.size()); 824 ShellLocation loc = sidesets[id - 1]; 825 element_surface_map(loc, elem_sides); 826 } 827 sideset_touching_blocks(int64_t)828 std::vector<std::string> GeneratedMesh::sideset_touching_blocks(int64_t /*set_id*/) const 829 { 830 std::vector<std::string> result(1, "block_1"); 831 return result; 832 } 833 set_variable_count(const std::string & type,size_t count)834 void GeneratedMesh::set_variable_count(const std::string &type, size_t count) 835 { 836 if (type == "global") { 837 variableCount[Ioss::REGION] = count; 838 } 839 else if (type == "element") { 840 variableCount[Ioss::ELEMENTBLOCK] = count; 841 } 842 else if (type == "nodal" || type == "node") { 843 variableCount[Ioss::NODEBLOCK] = count; 844 } 845 else if (type == "surface" || type == "sideset") { 846 variableCount[Ioss::SIDEBLOCK] = count; 847 } 848 else { 849 fmt::print(Ioss::WARNING(), 850 "(Iogs::GeneratedMesh::set_variable_count)\n" 851 " Unrecognized variable type '{}'. Valid types are:\n" 852 " global, element, node, nodal, surface, sideset.\n", 853 type); 854 } 855 } 856 set_rotation(const std::string & axis,double angle_degrees)857 void GeneratedMesh::set_rotation(const std::string &axis, double angle_degrees) 858 { 859 // PI / 180. Used in converting angle in degrees to radians 860 static double degang = std::atan2(0.0, -1.0) / 180.0; 861 862 doRotation = true; 863 864 int n1 = -1; 865 int n2 = -1; 866 int n3 = -1; 867 868 if (axis == "x" || axis == "X") { 869 n1 = 1; 870 n2 = 2; 871 n3 = 0; 872 } 873 else if (axis == "y" || axis == "Y") { 874 n1 = 2; 875 n2 = 0; 876 n3 = 1; 877 } 878 else if (axis == "z" || axis == "Z") { 879 n1 = 0; 880 n2 = 1; 881 n3 = 2; 882 } 883 else { 884 fmt::print("\nInvalid axis specification '{}'. Valid options are 'x', 'y', or 'z'\n", axis); 885 return; 886 } 887 888 double ang = angle_degrees * degang; // Convert angle in degrees to radians 889 double cosang = std::cos(ang); 890 double sinang = std::sin(ang); 891 892 assert(n1 >= 0 && n2 >= 0 && n3 >= 0); 893 std::array<std::array<double, 3>, 3> by; 894 by[n1][n1] = cosang; 895 by[n2][n1] = -sinang; 896 by[n1][n3] = 0.0; 897 by[n1][n2] = sinang; 898 by[n2][n2] = cosang; 899 by[n2][n3] = 0.0; 900 by[n3][n1] = 0.0; 901 by[n3][n2] = 0.0; 902 by[n3][n3] = 1.0; 903 904 std::array<std::array<double, 3>, 3> res; 905 for (int i = 0; i < 3; i++) { 906 res[i][0] = rotmat[i][0] * by[0][0] + rotmat[i][1] * by[1][0] + rotmat[i][2] * by[2][0]; 907 res[i][1] = rotmat[i][0] * by[0][1] + rotmat[i][1] * by[1][1] + rotmat[i][2] * by[2][1]; 908 res[i][2] = rotmat[i][0] * by[0][2] + rotmat[i][1] * by[1][2] + rotmat[i][2] * by[2][2]; 909 } 910 rotmat = res; 911 } 912 } // namespace Iogs 913