Sierra Toolkit  Version of the Day
GeneratedMesh.cpp
1 /*------------------------------------------------------------------------*/
2 /* Copyright 2010 Sandia Corporation. */
3 /* Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive */
4 /* license for use of this work by or on behalf of the U.S. Government. */
5 /* Export of this program may require a license from the */
6 /* United States Government. */
7 /*------------------------------------------------------------------------*/
8 
9 #include <stk_util/unit_test_support/GeneratedMesh.hpp>
10 #include <stk_util/util/tokenize.hpp>
11 
12 #include <cmath>
13 #include <cstring>
14 
15 #include <vector>
16 #include <cstdlib>
17 #include <string>
18 #include <iostream>
19 #include <iomanip>
20 #include <assert.h>
21 
22 namespace stk_classic {
23  namespace io {
24  namespace util {
25  GeneratedMesh::GeneratedMesh(int num_x, int num_y, int num_z,
26  int proc_count, int my_proc) :
27  numX(num_x), numY(num_y), numZ(num_z), myNumZ(num_z), myStartZ(0),
28  processorCount(proc_count), myProcessor(my_proc),
29  offX(0), offY(0), offZ(0),
30  sclX(1), sclY(1), sclZ(1),
31  doRotation(false)
32  {
33  initialize();
34  }
35 
36  GeneratedMesh::GeneratedMesh(const std::string &parameters,
37  int proc_count, int my_proc) :
38  numX(0), numY(0), numZ(0), myNumZ(0), myStartZ(0),
39  processorCount(proc_count), myProcessor(my_proc),
40  offX(0), offY(0), offZ(0),
41  sclX(1), sclY(1), sclZ(1),
42  doRotation(false)
43  {
44  std::vector<std::string> groups;
45  stk_classic::util::tokenize(parameters, "|+", groups);
46 
47  // First 'group' is the interval specification -- IxJxK
48  std::vector<std::string> tokens;
49  stk_classic::util::tokenize(groups[0], "x", tokens);
50  assert(tokens.size() == 3);
51  numX = std::strtol(tokens[0].c_str(), NULL, 10);
52  numY = std::strtol(tokens[1].c_str(), NULL, 10);
53  numZ = std::strtol(tokens[2].c_str(), NULL, 10);
54 
55  initialize();
56  parse_options(groups);
57 
58  }
59 
60  GeneratedMesh::~GeneratedMesh() {}
61 
62  void GeneratedMesh::initialize()
63  {
64  assert(numZ >= processorCount);
65  if (processorCount > 1) {
66  myNumZ = numZ / processorCount;
67  if (myProcessor < (numZ % processorCount)) myNumZ++;
68 
69  // Determine myStartZ for this processor...
70  size_t extra = numZ % processorCount;
71  if (extra > myProcessor)
72  extra = myProcessor;
73  size_t per_proc = numZ / processorCount;
74  myStartZ = myProcessor * per_proc + extra;
75  } else {
76  myNumZ = numZ;
77  }
78 
79  for (int i=0; i < 3; i++) {
80  for (int j=0; j < 3; j++) {
81  rotmat[i][j] = 0.0;
82  }
83  rotmat[i][i] = 1.0;
84  }
85  }
86 
87  size_t GeneratedMesh::add_shell_block(ShellLocation loc)
88  {
89  shellBlocks.push_back(loc);
90  return shellBlocks.size();
91  }
92 
93  size_t GeneratedMesh::add_nodeset(ShellLocation loc)
94  {
95  nodesets.push_back(loc);
96  return nodesets.size();
97  }
98 
99  size_t GeneratedMesh::add_sideset(ShellLocation loc)
100  {
101  sidesets.push_back(loc);
102  return sidesets.size();
103  }
104 
105  void GeneratedMesh::set_bbox(double xmin, double ymin, double zmin,
106  double xmax, double ymax, double zmax)
107  {
108  // NOTE: All calculations are based on the currently
109  // active interval settings. If scale or offset or zdecomp
110  // specified later in the option list, you may not get the
111  // desired bounding box.
112  double x_range = xmax - xmin;
113  double y_range = ymax - ymin;
114  double z_range = zmax - zmin;
115 
116  sclX = x_range / static_cast<double>(numX);
117  sclY = y_range / static_cast<double>(numY);
118  sclZ = z_range / static_cast<double>(numZ);
119 
120  offX = xmin;
121  offY = ymin;
122  offZ = zmin;
123  }
124 
125  void GeneratedMesh::set_scale(double scl_x, double scl_y, double scl_z)
126  {
127  sclX = scl_x;
128  sclY = scl_y;
129  sclZ = scl_z;
130  }
131 
132  void GeneratedMesh::set_offset(double off_x, double off_y, double off_z)
133  {
134  offX = off_x;
135  offY = off_y;
136  offZ = off_z;
137  }
138 
139  void GeneratedMesh::parse_options(const std::vector<std::string> &groups)
140  {
141  for (size_t i=1; i < groups.size(); i++) {
142  std::vector<std::string> option;
143  stk_classic::util::tokenize(groups[i], ":", option);
144  // option[0] is the type of the option and option[1] is the argument to the option.
145 
146  if (option[0] == "shell") {
147  // Option of the form "shell:xXyYzZ"
148  // The argument specifies whether there is a shell block
149  // at the location. 'x' is minX, 'X' is maxX, etc.
150  size_t length = option[1].size();
151  for (size_t j=0; j < length; j++) {
152  switch (option[1][j]) {
153  case 'x':
154  add_shell_block(MX);
155  break;
156  case 'X':
157  add_shell_block(PX);
158  break;
159  case 'y':
160  add_shell_block(MY);
161  break;
162  case 'Y':
163  add_shell_block(PY);
164  break;
165  case 'z':
166  add_shell_block(MZ);
167  break;
168  case 'Z':
169  add_shell_block(PZ);
170  break;
171  default:
172  std::cerr << "ERROR: Unrecognized shell location option '"
173  << option[1][j]
174  << "'.";
175  }
176  }
177  }
178  else if (option[0] == "nodeset") {
179  // Option of the form "nodeset:xXyYzZ"
180  // The argument specifies whether there is a nodeset
181  // at the location. 'x' is minX, 'X' is maxX, etc.
182  size_t length = option[1].size();
183  for (size_t j=0; j < length; j++) {
184  switch (option[1][j]) {
185  case 'x':
186  add_nodeset(MX);
187  break;
188  case 'X':
189  add_nodeset(PX);
190  break;
191  case 'y':
192  add_nodeset(MY);
193  break;
194  case 'Y':
195  add_nodeset(PY);
196  break;
197  case 'z':
198  add_nodeset(MZ);
199  break;
200  case 'Z':
201  add_nodeset(PZ);
202  break;
203  default:
204  std::cerr << "ERROR: Unrecognized nodeset location option '"
205  << option[1][j]
206  << "'.";
207  }
208  }
209  }
210  else if (option[0] == "sideset") {
211  // Option of the form "sideset:xXyYzZ"
212  // The argument specifies whether there is a sideset
213  // at the location. 'x' is minX, 'X' is maxX, etc.
214  size_t length = option[1].size();
215  for (size_t j=0; j < length; j++) {
216  switch (option[1][j]) {
217  case 'x':
218  add_sideset(MX);
219  break;
220  case 'X':
221  add_sideset(PX);
222  break;
223  case 'y':
224  add_sideset(MY);
225  break;
226  case 'Y':
227  add_sideset(PY);
228  break;
229  case 'z':
230  add_sideset(MZ);
231  break;
232  case 'Z':
233  add_sideset(PZ);
234  break;
235  default:
236  std::cerr << "ERROR: Unrecognized sideset location option '"
237  << option[1][j]
238  << "'.";
239  }
240  }
241  }
242  else if (option[0] == "scale") {
243  // Option of the form "scale:xs,ys,zs
244  std::vector<std::string> tokens;
245  stk_classic::util::tokenize(option[1], ",", tokens);
246  assert(tokens.size() == 3);
247  sclX = std::strtod(tokens[0].c_str(), NULL);
248  sclY = std::strtod(tokens[1].c_str(), NULL);
249  sclZ = std::strtod(tokens[2].c_str(), NULL);
250  }
251 
252  else if (option[0] == "offset") {
253  // Option of the form "offset:xo,yo,zo
254  std::vector<std::string> tokens;
255  stk_classic::util::tokenize(option[1], ",", tokens);
256  assert(tokens.size() == 3);
257  offX = std::strtod(tokens[0].c_str(), NULL);
258  offY = std::strtod(tokens[1].c_str(), NULL);
259  offZ = std::strtod(tokens[2].c_str(), NULL);
260  }
261 
262  else if (option[0] == "zdecomp") {
263  // Option of the form "zdecomp:1,1,2,2,1,2,...
264  // Specifies the number of intervals in the z direction
265  // for each processor. The number of tokens must match
266  // the number of processors. Note that the new numZ will
267  // be the sum of the intervals specified in this command.
268  std::vector<std::string> tokens;
269  stk_classic::util::tokenize(option[1], ",", tokens);
270  assert(tokens.size() == processorCount);
271  std::vector<size_t> Zs;
272  numZ = 0;
273  for (size_t j = 0; j < processorCount; j++) {
274  Zs.push_back(std::strtol(tokens[j].c_str(), NULL, 10));
275  numZ += Zs[j];
276  }
277  myNumZ = Zs[myProcessor];
278  myStartZ = 0;
279  for (size_t j=0; j < myProcessor; j++) {
280  myStartZ += Zs[j];
281  }
282  }
283 
284  else if (option[0] == "bbox") {
285  // Bounding-Box Option of the form "bbox:xmin,ymin,zmin,xmax,ymax,zmaxo
286  std::vector<std::string> tokens;
287  stk_classic::util::tokenize(option[1], ",", tokens);
288  assert(tokens.size() == 6);
289  double xmin = std::strtod(tokens[0].c_str(), NULL);
290  double ymin = std::strtod(tokens[1].c_str(), NULL);
291  double zmin = std::strtod(tokens[2].c_str(), NULL);
292  double xmax = std::strtod(tokens[3].c_str(), NULL);
293  double ymax = std::strtod(tokens[4].c_str(), NULL);
294  double zmax = std::strtod(tokens[5].c_str(), NULL);
295 
296  set_bbox(xmin, ymin, zmin, xmax, ymax, zmax);
297  }
298 
299  else if (option[0] == "rotate") {
300  // Rotate Option of the form "rotate:axis,angle,axis,angle,...
301  std::vector<std::string> tokens;
302  stk_classic::util::tokenize(option[1], ",", tokens);
303  assert(tokens.size() %2 == 0);
304  for (size_t ir=0; ir < tokens.size();) {
305  std::string axis = tokens[ir++];
306  double angle_degree = std::strtod(tokens[ir++].c_str(), NULL);
307  set_rotation(axis, angle_degree);
308  }
309  }
310 
311  else if (option[0] == "help") {
312  std::cerr << "\nValid Options for GeneratedMesh parameter string:\n"
313  << "\tIxJxK -- specifies intervals; must be first option. Ex: 4x10x12\n"
314  << "\toffset:xoff, yoff, zoff\n"
315  << "\tscale: xscl, yscl, zscl\n"
316  << "\tzdecomp:n1,n2,n3,...,n#proc\n"
317  << "\tbbox: xmin, ymin, zmin, xmax, ymax, zmax\n"
318  << "\trotate: axis,angle,axis,angle,...\n"
319  << "\tshell:xXyYzZ (specifies which plane to apply shell)\n"
320  << "\tnodeset:xXyXzZ (specifies which plane to apply nodeset)\n"
321  << "\tsideset:xXyXzZ (specifies which plane to apply sideset)\n"
322  << "\tshow -- show mesh parameters\n"
323  << "\thelp -- show this list\n\n";
324  }
325 
326  else if (option[0] == "show") {
327  show_parameters();
328  }
329 
330  else {
331  std::cerr << "ERROR: Unrecognized option '" << option[0]
332  << "'. It will be ignored.\n";
333  }
334  }
335  }
336 
337  void GeneratedMesh::show_parameters() const
338  {
339  if (myProcessor == 0) {
340  std::cerr << "\nMesh Parameters:\n"
341  << "\tIntervals: " << numX << " by " << numY << " by " << numZ << "\n"
342  << "\tX = " << sclX << " * (0.." << numX << ") + " << offX
343  << "\tRange: " << offX << " <= X <= " << offX + numX * sclX << "\n"
344  << "\tY = " << sclY << " * (0.." << numY << ") + " << offY
345  << "\tRange: " << offY << " <= Y <= " << offY + numY * sclY << "\n"
346  << "\tZ = " << sclZ << " * (0.." << numZ << ") + " << offZ
347  << "\tRange: " << offZ << " <= Z <= " << offZ + numZ * sclZ << "\n\n"
348  << "\tNode Count (total) = " << std::setw(9) << node_count() << "\n"
349  << "\tElement Count (total) = " << std::setw(9) << element_count() << "\n"
350  << "\tBlock Count = " << std::setw(9) << block_count() << "\n"
351  << "\tNodeset Count = " << std::setw(9) << nodeset_count() << "\n"
352  << "\tSideset Count = " << std::setw(9) << sideset_count() << "\n\n";
353  if (doRotation) {
354  std::cerr << "\tRotation Matrix: \n\t" << std::scientific ;
355  for (int ii=0; ii < 3; ii++) {
356  for (int jj=0; jj < 3; jj++) {
357  std::cerr << std::setw(14) << rotmat[ii][jj] << "\t";
358  }
359  std::cerr << "\n\t";
360  }
361  std::cerr << std::fixed << "\n";
362  }
363  }
364  }
365 
366  size_t GeneratedMesh::node_count() const
367  {
368  return (numX+1) * (numY+1) * (numZ+1);
369  }
370 
371  size_t GeneratedMesh::node_count_proc() const
372  {
373  return (numX+1) * (numY+1) * (myNumZ+1);
374  }
375 
376  size_t GeneratedMesh::block_count() const
377  {
378  return shellBlocks.size() + 1;
379  }
380 
381  size_t GeneratedMesh::nodeset_count() const
382  {
383  return nodesets.size();
384  }
385 
386  size_t GeneratedMesh::sideset_count() const
387  {
388  return sidesets.size();
389  }
390 
391  size_t GeneratedMesh::element_count() const
392  {
393  size_t count = element_count(1);
394  for (size_t i=0; i < shellBlocks.size(); i++) {
395  count += element_count(i+2);
396  }
397  return count;
398  }
399 
400  size_t GeneratedMesh::element_count_proc() const
401  {
402  size_t count = 0;
403  for (size_t i=0; i < block_count(); i++) {
404  count += element_count_proc(i+1);
405  }
406  return count;
407  }
408 
409  size_t GeneratedMesh::element_count(size_t block_number) const
410  {
411  assert(block_number <= block_count());
412 
413  if (block_number == 1) {
414  return numX * numY * numZ;
415  } else {
416  ShellLocation loc = shellBlocks[block_number-2];
417  return shell_element_count(loc);
418  }
419  }
420 
421  size_t GeneratedMesh::shell_element_count(ShellLocation loc) const
422  {
423  switch (loc) {
424  case MX:
425  case PX:
426  return numY * numZ;
427  case MY:
428  case PY:
429  return numX * numZ;
430  case MZ:
431  case PZ:
432  return numX * numY;
433  }
434  return 0;
435  }
436 
437  size_t GeneratedMesh::element_count_proc(size_t block_number) const
438  {
439  assert(block_number <= block_count());
440 
441  if (block_number == 1) {
442  return numX * numY * myNumZ;
443  } else {
444  ShellLocation loc = shellBlocks[block_number-2];
445  return shell_element_count_proc(loc);
446  }
447  }
448 
449  size_t GeneratedMesh::shell_element_count_proc(ShellLocation loc) const
450  {
451  switch (loc) {
452  case MX:
453  case PX:
454  return numY * myNumZ;
455  case MY:
456  case PY:
457  return numX * myNumZ;
458  case MZ:
459  if (myProcessor == 0)
460  return numX * numY;
461  else
462  return 0;
463  case PZ:
464  if (myProcessor == processorCount -1)
465  return numX * numY;
466  else
467  return 0;
468  }
469  return 0;
470  }
471 
472  size_t GeneratedMesh::nodeset_node_count(size_t id) const
473  {
474  // id is position in nodeset list + 1
475  assert(id > 0 && id <= nodesets.size());
476  ShellLocation loc = nodesets[id-1];
477  switch (loc) {
478  case MX:
479  case PX:
480  return (numY+1) * (numZ+1);
481  case MY:
482  case PY:
483  return (numX+1) * (numZ+1);
484  case MZ:
485  case PZ:
486  return (numX+1) * (numY+1);
487  }
488  return 0;
489  }
490 
491  size_t GeneratedMesh::nodeset_node_count_proc(size_t id) const
492  {
493  // id is position in nodeset list + 1
494  assert(id > 0 && id <= nodesets.size());
495  ShellLocation loc = nodesets[id-1];
496  switch (loc) {
497  case MX:
498  case PX:
499  return (numY+1) * (myNumZ+1);
500  case MY:
501  case PY:
502  return (numX+1) * (myNumZ+1);
503  case MZ:
504  if (myProcessor == 0)
505  return (numX+1) * (numY+1);
506  else
507  return 0;
508  case PZ:
509  if (myProcessor == processorCount -1)
510  return (numX+1) * (numY+1);
511  else
512  return 0;
513  }
514  return 0;
515  }
516 
517  size_t GeneratedMesh::sideset_side_count(size_t id) const
518  {
519  // id is position in sideset list + 1
520  assert(id > 0 && id <= sidesets.size());
521  ShellLocation loc = sidesets[id-1];
522  switch (loc) {
523  case MX:
524  case PX:
525  return numY * numZ;
526  case MY:
527  case PY:
528  return numX * numZ;
529  case MZ:
530  case PZ:
531  return numX * numY;
532  }
533  return 0;
534  }
535 
536  size_t GeneratedMesh::sideset_side_count_proc(size_t id) const
537  {
538  // id is position in sideset list + 1
539  assert(id > 0 && id <= sidesets.size());
540  ShellLocation loc = sidesets[id-1];
541  switch (loc) {
542  case MX:
543  case PX:
544  return numY * myNumZ;
545  case MY:
546  case PY:
547  return numX * myNumZ;
548  case MZ:
549  if (myProcessor == 0)
550  return numX * numY;
551  else
552  return 0;
553  case PZ:
554  if (myProcessor == processorCount -1)
555  return numX * numY;
556  else
557  return 0;
558  }
559  return 0;
560  }
561 
562  std::pair<std::string,int> GeneratedMesh::topology_type(size_t block_number) const
563  {
564  assert(block_number <= block_count() && block_number > 0);
565 
566  if (block_number == 1) {
567  return std::make_pair(std::string("hex8"), 8);
568  } else {
569  return std::make_pair(std::string("shell4"), 4);
570  }
571  }
572 
573  void GeneratedMesh::node_map(std::vector<int> &map)
574  {
575  size_t count = node_count_proc();
576  map.resize(count);
577  size_t offset = myStartZ * (numX+1) * (numY+1);
578  for (size_t i=0; i < count; i++) {
579  map[i] = static_cast<int>(offset + i + 1);
580  }
581  }
582 
583  size_t GeneratedMesh::communication_node_count_proc() const
584  {
585  size_t count = (numX+1) * (numY+1);
586  if (myProcessor != 0 && myProcessor != processorCount-1)
587  count *= 2;
588 
589  return count;
590  }
591 
592  void GeneratedMesh::node_communication_map(std::vector<int> &map, std::vector<int> &proc)
593  {
594  size_t count = (numX+1) * (numY+1);
595  size_t slab = count;
596  if (myProcessor != 0 && myProcessor != processorCount-1)
597  count *= 2;
598 
599  map.resize(count);
600  proc.resize(count);
601  size_t j = 0;
602  if (myProcessor != 0) {
603  size_t offset = myStartZ * (numX+1) * (numY+1);
604  for (size_t i=0; i < slab; i++) {
605  map[j] = static_cast<int>(offset + i + 1);
606  proc[j++] = static_cast<int>(myProcessor-1);
607  }
608  }
609  if (myProcessor != processorCount-1) {
610  size_t offset = (myStartZ + myNumZ) * (numX+1) * (numY+1);
611  for (size_t i=0; i < slab; i++) {
612  map[j] = static_cast<int>(offset + i + 1);
613  proc[j++] = static_cast<int>(myProcessor+1);
614  }
615  }
616  }
617 
618  void GeneratedMesh::element_map(size_t block_number, std::vector<int> &map) const
619  {
620  assert(block_number <= block_count() && block_number > 0);
621 
622  size_t count = element_count_proc(block_number);
623  map.resize(count);
624 
625  if (block_number == 1) {
626  // Hex block...
627  count = element_count_proc(1);
628  size_t offset = myStartZ * numX * numY;
629  for (size_t i=0; i < count; i++) {
630  map[i] = static_cast<int>(offset + i + 1);
631  }
632  } else {
633  size_t start = element_count(1);
634 
635  // Shell blocks...
636  for (size_t ib=0; ib < shellBlocks.size(); ib++) {
637  count = element_count_proc(ib+2);
638  if (static_cast<size_t>(block_number) == ib + 2) {
639  size_t offset = 0;
640  ShellLocation loc = shellBlocks[ib];
641 
642  switch (loc) {
643  case MX:
644  case PX:
645  offset = myStartZ * numY;
646  break;
647 
648  case MY:
649  case PY:
650  offset = myStartZ * numX;
651  break;
652 
653  case MZ:
654  case PZ:
655  offset = 0;
656  break;
657  }
658  for (size_t i=0; i < count; i++) {
659  map[i] = static_cast<int>(start + offset + i + 1);
660  }
661  } else {
662  start += element_count(ib+2);
663  }
664  }
665  }
666  }
667 
668  void GeneratedMesh::element_map(std::vector<int> &map) const
669  {
670  size_t count = element_count_proc();
671  map.resize(count);
672 
673  size_t k = 0;
674  // Hex block...
675  count = element_count_proc(1);
676  size_t offset = myStartZ * numX * numY;
677  for (size_t i=0; i < count; i++) {
678  map[k++] = static_cast<int>(offset + i + 1);
679  }
680 
681  size_t start = element_count(1);
682 
683  // Shell blocks...
684  for (size_t ib=0; ib < shellBlocks.size(); ib++) {
685  count = element_count_proc(ib+2);
686  offset = 0;
687  ShellLocation loc = shellBlocks[ib];
688 
689  switch (loc) {
690  case MX:
691  case PX:
692  offset = myStartZ * numY;
693  break;
694 
695  case MY:
696  case PY:
697  offset = myStartZ * numX;
698  break;
699 
700  case MZ:
701  case PZ:
702  offset = 0;
703  break;
704  }
705  for (size_t i=0; i < count; i++) {
706  map[k++] = static_cast<int>(start + offset + i + 1);
707  }
708  start += element_count(ib+2);
709  }
710  }
711 
712  void GeneratedMesh::element_surface_map(ShellLocation loc, std::vector<int> &map) const
713  {
714  size_t count = shell_element_count_proc(loc);
715  map.resize(2*count);
716 
717  size_t index = 0;
718  size_t offset = 0;
719 
720  switch (loc) {
721  case MX:
722  offset = myStartZ * numX * numY + 1; // 1-based elem id
723  for (size_t k = 0; k < myNumZ; ++k) {
724  for (size_t j = 0; j < numY; ++j) {
725  map[index++] = offset;
726  map[index++] = 3; // 0-based local face id
727  offset += numX;
728  }
729  }
730  break;
731 
732  case PX:
733  offset = myStartZ * numX * numY + numX;
734  for (size_t k = 0; k < myNumZ; ++k) {
735  for (size_t j = 0; j < numY; ++j) {
736  map[index++] = offset; // 1-based elem id
737  map[index++] = 1; // 0-based local face id
738  offset += numX;
739  }
740  }
741  break;
742 
743  case MY:
744  offset = myStartZ * numX * numY + 1;
745  for (size_t k = 0; k < myNumZ; ++k) {
746  for (size_t i = 0; i < numX; ++i) {
747  map[index++] = offset++;
748  map[index++] = 0; // 0-based local face id
749  }
750  offset+= numX * (numY-1);
751  }
752  break;
753 
754  case PY:
755  offset = myStartZ * numX * numY + numX * (numY-1) +1;
756  for (size_t k = 0; k < myNumZ; ++k) {
757  for (size_t i = 0; i < numX; ++i) {
758  map[index++] = offset++;
759  map[index++] = 2; // 0-based local face id
760  }
761  offset+= numX * (numY-1);
762  }
763  break;
764 
765  case MZ:
766  if (myProcessor == 0) {
767  offset = 1;
768  for (size_t i=0; i < numY; i++) {
769  for (size_t j=0; j < numX; j++) {
770  map[index++] = offset++;
771  map[index++] = 4;
772  }
773  }
774  }
775  break;
776 
777  case PZ:
778  if (myProcessor == processorCount-1) {
779  offset = (numZ-1)*numX*numY + 1;
780  for (size_t i=0, k=0; i < numY; i++) {
781  for (size_t j=0; j < numX; j++, k++) {
782  map[index++] = offset++;
783  map[index++] = 5;
784  }
785  }
786  }
787  break;
788  }
789  }
790 
791  void GeneratedMesh::coordinates(std::vector<double> &coord) const
792  {
793  /* create global coordinates */
794  size_t count = node_count_proc();
795  coord.resize(count * 3);
796 
797  size_t k = 0;
798  for (size_t m=myStartZ; m < myStartZ+myNumZ+1; m++) {
799  for (size_t i=0; i < numY+1; i++) {
800  for (size_t j=0; j < numX+1; j++) {
801  coord[k++] = sclX * static_cast<double>(j) + offX;
802  coord[k++] = sclY * static_cast<double>(i) + offY;
803  coord[k++] = sclZ * static_cast<double>(m) + offZ;
804  }
805  }
806  }
807 
808  if (doRotation) {
809  for (size_t i=0; i < count*3; i+=3) {
810  double xn = coord[i+0];
811  double yn = coord[i+1];
812  double zn = coord[i+2];
813  coord[i+0] = xn * rotmat[0][0] + yn * rotmat[1][0] + zn * rotmat[2][0];
814  coord[i+1] = xn * rotmat[0][1] + yn * rotmat[1][1] + zn * rotmat[2][1];
815  coord[i+2] = xn * rotmat[0][2] + yn * rotmat[1][2] + zn * rotmat[2][2];
816  }
817  }
818  }
819 
820  void GeneratedMesh::coordinates(std::vector<double> &x,
821  std::vector<double> &y,
822  std::vector<double> &z) const
823  {
824  /* create global coordinates */
825  size_t count = node_count_proc();
826  x.resize(count);
827  y.resize(count);
828  z.resize(count);
829 
830  size_t k = 0;
831  for (size_t m=myStartZ; m < myStartZ+myNumZ+1; m++) {
832  for (size_t i=0; i < numY+1; i++) {
833  for (size_t j=0; j < numX+1; j++) {
834  x[k] = sclX * static_cast<double>(j) + offX;
835  y[k] = sclY * static_cast<double>(i) + offY;
836  z[k] = sclZ * static_cast<double>(m) + offZ;
837  ++k;
838  }
839  }
840  }
841  if (doRotation) {
842  for (size_t i=0; i < count; i++) {
843  double xn = x[i];
844  double yn = y[i];
845  double zn = z[i];
846  x[i] = xn * rotmat[0][0] + yn * rotmat[1][0] + zn * rotmat[2][0];
847  y[i] = xn * rotmat[0][1] + yn * rotmat[1][1] + zn * rotmat[2][1];
848  z[i] = xn * rotmat[0][2] + yn * rotmat[1][2] + zn * rotmat[2][2];
849  }
850  }
851  }
852 
853  void GeneratedMesh::connectivity(size_t block_number, std::vector<int> &connect) const
854  {
855  assert(block_number <= block_count());
856 
857  size_t xp1yp1 = (numX+1) * (numY+1);
858 
859  /* build connectivity array (node list) for mesh */
860  if (block_number == 1) { // HEX Element Block
861  connect.resize(element_count_proc(block_number)*8);
862 
863  size_t cnt = 0;
864  for (size_t m=myStartZ; m < myNumZ+myStartZ; m++) {
865  for (size_t i=0, k=0; i < numY; i++) {
866  for (size_t j=0; j < numX; j++, k++) {
867  size_t base = (m*xp1yp1) + k + i + 1;
868  ;
869  connect[cnt++] = static_cast<int>(base);
870  connect[cnt++] = static_cast<int>(base+1);
871  connect[cnt++] = static_cast<int>(base+numX+2);
872  connect[cnt++] = static_cast<int>(base+numX+1);
873 
874  connect[cnt++] = static_cast<int>(xp1yp1 + base);
875  connect[cnt++] = static_cast<int>(xp1yp1 + base+1);
876  connect[cnt++] = static_cast<int>(xp1yp1 + base+numX+2);
877  connect[cnt++] = static_cast<int>(xp1yp1 + base+numX+1);
878  }
879  }
880  }
881  } else { // Shell blocks....
882  ShellLocation loc = shellBlocks[block_number-2];
883  connect.resize(element_count_proc(block_number)*4);
884 
885  size_t cnt = 0;
886  switch (loc) {
887  case MX: // Minumum X Face
888  for (size_t i=0; i < myNumZ; i++) {
889  size_t layer_off = i * xp1yp1;
890  for (size_t j=0; j < numY; j++) {
891  size_t base = layer_off + j * (numX+1) + 1 + myStartZ * xp1yp1;
892  connect[cnt++] = static_cast<int>(base);
893  connect[cnt++] = static_cast<int>(base + xp1yp1);
894  connect[cnt++] = static_cast<int>(base + xp1yp1 + (numX+1));
895  connect[cnt++] = static_cast<int>(base + (numX+1));
896  }
897  }
898  break;
899  case PX: // Maximum X Face
900  for (size_t i=0; i < myNumZ; i++) {
901  size_t layer_off = i * xp1yp1;
902  for (size_t j=0; j < numY; j++) {
903  size_t base = layer_off + j * (numX+1) + numX + 1 + myStartZ * xp1yp1;
904  connect[cnt++] = static_cast<int>(base);
905  connect[cnt++] = static_cast<int>(base + (numX+1));
906  connect[cnt++] = static_cast<int>(base + xp1yp1 + (numX+1));
907  connect[cnt++] = static_cast<int>(base + xp1yp1);
908  }
909  }
910  break;
911  case MY: // Minumum Y Face
912  for (size_t i=0; i < myNumZ; i++) {
913  size_t layer_off = i * xp1yp1;
914  for (size_t j=0; j < numX; j++) {
915  size_t base = layer_off + j + 1 + myStartZ * xp1yp1;
916  connect[cnt++] = static_cast<int>(base);
917  connect[cnt++] = static_cast<int>(base + 1);
918  connect[cnt++] = static_cast<int>(base + xp1yp1 + 1);
919  connect[cnt++] = static_cast<int>(base + xp1yp1);
920  }
921  }
922  break;
923  case PY: // Maximum Y Face
924  for (size_t i=0; i < myNumZ; i++) {
925  size_t layer_off = i * xp1yp1;
926  for (size_t j=0; j < numX; j++) {
927  size_t base = layer_off + (numX+1)*(numY) + j + 1 + myStartZ * xp1yp1;
928  connect[cnt++] = static_cast<int>(base);
929  connect[cnt++] = static_cast<int>(base + xp1yp1);
930  connect[cnt++] = static_cast<int>(base + xp1yp1 + 1);
931  connect[cnt++] = static_cast<int>(base + 1);
932  }
933  }
934  break;
935  case MZ: // Minumum Z Face
936  if (myProcessor == 0) {
937  for (size_t i=0, k=0; i < numY; i++) {
938  for (size_t j=0; j < numX; j++, k++) {
939  size_t base = i + k + 1 + myStartZ * xp1yp1;
940  connect[cnt++] = static_cast<int>(base);
941  connect[cnt++] = static_cast<int>(base+numX+1);
942  connect[cnt++] = static_cast<int>(base+numX+2);
943  connect[cnt++] = static_cast<int>(base+1);
944  }
945  }
946  }
947  break;
948  case PZ: // Maximum Z Face
949  if (myProcessor == processorCount-1) {
950  for (size_t i=0, k=0; i < numY; i++) {
951  for (size_t j=0; j < numX; j++, k++) {
952  size_t base = xp1yp1 * (numZ - myStartZ) + k + i + 1 + myStartZ * xp1yp1;
953  connect[cnt++] = static_cast<int>(base);
954  connect[cnt++] = static_cast<int>(base+1);
955  connect[cnt++] = static_cast<int>(base+numX+2);
956  connect[cnt++] = static_cast<int>(base+numX+1);
957  }
958  }
959  }
960  break;
961  }
962  assert(cnt == 4 * element_count_proc(block_number));
963  }
964  return;
965  }
966 
967  void GeneratedMesh::nodeset_nodes(size_t id, std::vector<int> &nodes) const
968  {
969  // id is position in nodeset list + 1
970  assert(id > 0 && id <= nodesets.size());
971  ShellLocation loc = nodesets[id-1];
972  nodes.resize(nodeset_node_count_proc(id));
973 
974  size_t xp1yp1 = (numX+1) * (numY+1);
975  size_t k = 0;
976 
977  switch (loc) {
978  case MX: // Minumum X Face
979  for (size_t i=0; i < myNumZ+1; i++) {
980  size_t layer_off = myStartZ * xp1yp1 + i * xp1yp1;
981  for (size_t j=0; j < numY+1; j++) {
982  nodes[k++] = layer_off + j * (numX+1) + 1;
983  }
984  }
985  break;
986  case PX: // Maximum X Face
987  for (size_t i=0; i < myNumZ+1; i++) {
988  size_t layer_off = myStartZ * xp1yp1 + i * xp1yp1;
989  for (size_t j=0; j < numY+1; j++) {
990  nodes[k++] = layer_off + j * (numX+1) + numX + 1;
991  }
992  }
993  break;
994  case MY: // Minumum Y Face
995  for (size_t i=0; i < myNumZ+1; i++) {
996  size_t layer_off = myStartZ * xp1yp1 + i * xp1yp1;
997  for (size_t j=0; j < numX+1; j++) {
998  nodes[k++] = layer_off + j + 1;
999  }
1000  }
1001  break;
1002  case PY: // Maximum Y Face
1003  for (size_t i=0; i < myNumZ+1; i++) {
1004  size_t layer_off = myStartZ * xp1yp1 + i * xp1yp1;
1005  for (size_t j=0; j < numX+1; j++) {
1006  nodes[k++] = layer_off + (numX+1)*(numY) + j + 1;
1007  }
1008  }
1009  break;
1010  case MZ: // Minumum Z Face
1011  if (myProcessor == 0) {
1012  for (size_t i=0; i < (numY+1) * (numX+1); i++) {
1013  nodes[i] = i+1;
1014  }
1015  }
1016  break;
1017  case PZ: // Maximum Z Face
1018  if (myProcessor == processorCount-1) {
1019  size_t offset = (numY+1) * (numX+1) * numZ;
1020  for (size_t i=0; i < (numY+1) * (numX+1); i++) {
1021  nodes[i] = offset + i+1;
1022  }
1023  }
1024  break;
1025  }
1026  }
1027 
1028  void GeneratedMesh::sideset_elem_sides(size_t id, std::vector<int> &elem_sides) const
1029  {
1030  // id is position in sideset list + 1
1031  assert(id > 0 && id <= sidesets.size());
1032  ShellLocation loc = sidesets[id-1];
1033 
1034  // If there is a shell block on this face, then the sideset is
1035  // applied to the shell block; if not, it is applied to the
1036  // underlying hex elements.
1037  bool underlying_shell = false;
1038  size_t shell_block = 0;
1039  for (size_t i = 0; i < shellBlocks.size(); i++) {
1040  if (shellBlocks[i] == loc) {
1041  underlying_shell = true;
1042  shell_block = i+2;
1043  break;
1044  }
1045  }
1046 
1047  if (underlying_shell) {
1048  // Get ids of shell elements at this location...
1049  element_map(shell_block, elem_sides);
1050 
1051  // Insert face_ordinal in between each entry in elem_sides...
1052  // Face will be 0 for all shells...
1053  elem_sides.resize(2*sideset_side_count_proc(id));
1054  int face_ordinal = 0;
1055  int i = 2* (int)sideset_side_count_proc(id) - 1;
1056  int j = (int)sideset_side_count_proc(id) - 1;
1057  while (i >= 0) {
1058  elem_sides[i--] = face_ordinal;
1059  elem_sides[i--] = elem_sides[j--];
1060  }
1061  } else {
1062  element_surface_map(loc, elem_sides);
1063  }
1064  }
1065 
1066  void GeneratedMesh::set_rotation(const std::string &axis, double angle_degrees)
1067  {
1068  // PI / 180. Used in converting angle in degrees to radians
1069  static double degang = std::atan2(0.0, -1.0) / 180.0;
1070 
1071  doRotation = true;
1072 
1073  int n1 = -1;
1074  int n2 = -1;
1075  int n3 = -1;
1076 
1077  if (axis == "x" || axis == "X") {
1078  n1 = 1; n2 = 2; n3 = 0;
1079  } else if (axis == "y" || axis == "Y") {
1080  n1 = 2; n2 = 0; n3 = 1;
1081  } else if (axis == "z" || axis == "Z") {
1082  n1 = 0; n2 = 1; n3 = 2;
1083  } else {
1084  std::cerr << "\nInvalid axis specification '" << axis << "'. Valid options are 'x', 'y', or 'z'\n";
1085  return;
1086  }
1087 
1088  double ang = angle_degrees * degang; // Convert angle in degrees to radians
1089  double cosang = std::cos(ang);
1090  double sinang = std::sin(ang);
1091 
1092  assert(n1 >= 0 && n2 >= 0 && n3 >= 0);
1093  double by[3][3];
1094  by[n1][n1] = cosang;
1095  by[n2][n1] = -sinang;
1096  by[n1][n3] = 0.0;
1097  by[n1][n2] = sinang;
1098  by[n2][n2] = cosang;
1099  by[n2][n3] = 0.0;
1100  by[n3][n1] = 0.0;
1101  by[n3][n2] = 0.0;
1102  by[n3][n3] = 1.0;
1103 
1104  double res[3][3];
1105  for (size_t i=0; i < 3; i++) {
1106  res[i][0] = rotmat[i][0]*by[0][0] + rotmat[i][1]*by[1][0] + rotmat[i][2]*by[2][0];
1107  res[i][1] = rotmat[i][0]*by[0][1] + rotmat[i][1]*by[1][1] + rotmat[i][2]*by[2][1];
1108  res[i][2] = rotmat[i][0]*by[0][2] + rotmat[i][1]*by[1][2] + rotmat[i][2]*by[2][2];
1109  }
1110 
1111 #if 1
1112  std::memcpy(rotmat, res, 9*sizeof(double));
1113 #else
1114  for (int i=0; i < 3; i++) {
1115  for (int j=0; j < 3; j++) {
1116  rotmat[i][j] = res[i][j];
1117  }
1118  }
1119 #endif
1120  }
1121  }
1122  }
1123 }
1124 
1125 #if defined(DEBUG)
1126 #include <sstream>
1127  std::string ToString(size_t t) {
1128  std::ostringstream os;
1129  os << t;
1130  return os.str();
1131  }
1132 
1133 #include <exodusII.h>
1134 
1135  int main() {
1136  int num_processors = 8;
1137  for (int proc = 0; proc < num_processors; proc++) {
1138 
1139  stk_classic::GeneratedMesh mesh(100, 125, 10*num_processors, num_processors, proc);
1140 
1141  std::cerr << "Node Count (total) = " << mesh.node_count() << "\n";
1142  std::cerr << "Node Count (proc) = " << mesh.node_count_proc() << "\n";
1143  std::cerr << "Element Count (total) = " << mesh.element_count() << "\n";
1144  std::cerr << "Element Count (proc) = " << mesh.element_count_proc() << "\n";
1145  std::cerr << "Block Count = " << mesh.block_count() << "\n";
1146 
1147  int CPU_word_size = 8; /* sizeof(float) */
1148  int IO_word_size = 8; /* (4 bytes) */
1149  std::string name = "test-scale.e";
1150  if (num_processors > 1) {
1151  name += "." + ToString(num_processors) + "." + ToString(proc);
1152  }
1153  int exoid = ex_create (name.c_str(), EX_CLOBBER, &CPU_word_size, &IO_word_size);
1154 
1155  int num_nodes = mesh.node_count_proc();
1156  int num_elems = mesh.element_count_proc();
1157  int num_elem_blk = mesh.block_count();
1158  int error = ex_put_init (exoid, "title", 3, num_nodes, num_elems,
1159  num_elem_blk, 0, 0);
1160 
1161  if (num_processors > 1) {
1162  std::vector<int> nodes;
1163  std::vector<int> procs;
1164  mesh.node_communication_map(nodes, procs);
1165 
1166  int node_map_ids[1] = {1};
1167  int node_map_node_cnts[1] = {procs.size()};
1168  ex_put_init_info(exoid, num_processors, 1, "p");
1169  ex_put_loadbal_param(exoid, 0, 0, 0, 0, 0, 1, 0, proc);
1170  ex_put_cmap_params(exoid, node_map_ids, node_map_node_cnts, 0, 0, proc);
1171  ex_put_node_cmap(exoid, 1, &nodes[0], &procs[0], proc);
1172  }
1173 
1174  for (int i=1; i < mesh.block_count(); i++) {
1175  std::cerr << "Block " << i+1 << " has " << mesh.element_count_proc(i+1) << " "
1176  << mesh.topology_type(i+1).first <<" elements\n";
1177 
1178  std::string btype = mesh.topology_type(i+1).first;
1179  error = ex_put_elem_block(exoid, i+1, btype.c_str(),
1180  mesh.element_count_proc(i+1),
1181  mesh.topology_type(i+1).second, 0);
1182  }
1183  {
1184  std::cerr << "Block " << 1 << " has " << mesh.element_count_proc(1) << " "
1185  << mesh.topology_type(1).first <<" elements\n";
1186 
1187  std::string btype = mesh.topology_type(1).first;
1188  error = ex_put_elem_block(exoid, 1, btype.c_str(),
1189  mesh.element_count_proc(1),
1190  mesh.topology_type(1).second, 0);
1191  }
1192 
1193  if (num_processors > 1) {
1194  std::vector<int> map;
1195  mesh.node_map(map);
1196  ex_put_id_map(exoid, EX_NODE_MAP, &map[0]);
1197 
1198  mesh.element_map(map);
1199  ex_put_id_map(exoid, EX_ELEM_MAP, &map[0]);
1200  }
1201 
1202  std::cerr << "Outputting connectivity...\n";
1203  for (int i=1; i < mesh.block_count(); i++) {
1204  if (mesh.element_count_proc(i+1) > 0) {
1205  std::vector<int> connectivity;
1206  mesh.connectivity(i+1, connectivity);
1207  ex_put_elem_conn(exoid, i+1, &connectivity[0]);
1208  }
1209  }
1210  {
1211  std::vector<int> connectivity;
1212  mesh.connectivity(1, connectivity);
1213  ex_put_elem_conn(exoid, 1, &connectivity[0]);
1214  }
1215 
1216  {
1217  std::vector<double> x;
1218  std::vector<double> y;
1219  std::vector<double> z;
1220  mesh.coordinates(x, y, z);
1221  error = ex_put_coord (exoid, &x[0], &y[0], &z[0]);
1222  }
1223 
1224  ex_close(exoid);
1225  }
1226  }
1227 #endif
_setw setw(int width)
Function setw sets the width for the next field as a manipulator.
Definition: WriterManip.hpp:44
Sierra Toolkit.
eastl::iterator_traits< InputIterator >::difference_type count(InputIterator first, InputIterator last, const T &value)