Sierra Toolkit Version of the Day
GeneratedMesh.cpp
00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010 Sandia Corporation.                     */
00003 /*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
00004 /*  license for use of this work by or on behalf of the U.S. Government.  */
00005 /*  Export of this program may require a license from the                 */
00006 /*  United States Government.                                             */
00007 /*------------------------------------------------------------------------*/
00008 
00009 #include <stk_util/unit_test_support/GeneratedMesh.hpp>
00010 #include <stk_util/util/tokenize.hpp>
00011 
00012 #include <cmath>
00013 #include <cstring>
00014 
00015 #include <vector>
00016 #include <cstdlib>
00017 #include <string>
00018 #include <iostream>
00019 #include <iomanip>
00020 #include <assert.h>
00021 
00022 namespace stk {
00023   namespace io {
00024     namespace util {
00025       GeneratedMesh::GeneratedMesh(int num_x, int num_y, int num_z,
00026                                    int proc_count, int my_proc) :
00027         numX(num_x), numY(num_y), numZ(num_z), myNumZ(num_z), myStartZ(0),
00028         processorCount(proc_count), myProcessor(my_proc),
00029         offX(0), offY(0), offZ(0),
00030         sclX(1), sclY(1), sclZ(1),
00031         doRotation(false)
00032       {
00033         initialize();
00034       }
00035 
00036       GeneratedMesh::GeneratedMesh(const std::string &parameters,
00037                                    int proc_count, int my_proc) :
00038         numX(0), numY(0), numZ(0), myNumZ(0), myStartZ(0),
00039         processorCount(proc_count), myProcessor(my_proc),
00040         offX(0), offY(0), offZ(0),
00041         sclX(1), sclY(1), sclZ(1),
00042         doRotation(false)
00043       {
00044         std::vector<std::string> groups;
00045         stk::util::tokenize(parameters, "|+", groups);
00046 
00047         // First 'group' is the interval specification -- IxJxK
00048         std::vector<std::string> tokens;
00049         stk::util::tokenize(groups[0], "x", tokens);
00050         assert(tokens.size() == 3);
00051         numX = std::strtol(tokens[0].c_str(), NULL, 10);
00052         numY = std::strtol(tokens[1].c_str(), NULL, 10);
00053         numZ = std::strtol(tokens[2].c_str(), NULL, 10);
00054 
00055         initialize();
00056         parse_options(groups);
00057 
00058       }
00059 
00060       GeneratedMesh::~GeneratedMesh() {}
00061 
00062       void GeneratedMesh::initialize()
00063       {
00064         assert(numZ >= processorCount);
00065         if (processorCount > 1) {
00066           myNumZ = numZ / processorCount;
00067           if (myProcessor < (numZ % processorCount)) myNumZ++;
00068 
00069           // Determine myStartZ for this processor...
00070           size_t extra = numZ % processorCount;
00071           if (extra > myProcessor)
00072             extra = myProcessor;
00073           size_t per_proc  = numZ / processorCount;
00074           myStartZ = myProcessor * per_proc + extra;
00075         } else {
00076           myNumZ = numZ;
00077         }
00078 
00079         for (int i=0; i < 3; i++) {
00080           for (int j=0; j < 3; j++) {
00081             rotmat[i][j] = 0.0;
00082           }
00083           rotmat[i][i] = 1.0;
00084         }
00085       }
00086 
00087       size_t GeneratedMesh::add_shell_block(ShellLocation loc)
00088       {
00089         shellBlocks.push_back(loc);
00090         return shellBlocks.size();
00091       }
00092 
00093       size_t GeneratedMesh::add_nodeset(ShellLocation loc)
00094       {
00095   nodesets.push_back(loc);
00096   return nodesets.size();
00097       }
00098       
00099       size_t GeneratedMesh::add_sideset(ShellLocation loc)
00100       {
00101   sidesets.push_back(loc);
00102   return sidesets.size();
00103       }
00104 
00105       void GeneratedMesh::set_bbox(double xmin, double ymin, double zmin,
00106                                    double xmax, double ymax, double zmax)
00107       {
00108         // NOTE: All calculations are based on the currently
00109         // active interval settings. If scale or offset or zdecomp
00110         // specified later in the option list, you may not get the
00111         // desired bounding box.
00112         double x_range = xmax - xmin;
00113         double y_range = ymax - ymin;
00114         double z_range = zmax - zmin;
00115 
00116         sclX = x_range / static_cast<double>(numX);
00117         sclY = y_range / static_cast<double>(numY);
00118         sclZ = z_range / static_cast<double>(numZ);
00119 
00120         offX = xmin;
00121         offY = ymin;
00122         offZ = zmin;
00123       }
00124 
00125       void GeneratedMesh::set_scale(double scl_x, double scl_y, double scl_z)
00126       {
00127         sclX = scl_x;
00128         sclY = scl_y;
00129         sclZ = scl_z;
00130       }
00131 
00132       void GeneratedMesh::set_offset(double off_x, double off_y, double off_z)
00133       {
00134         offX = off_x;
00135         offY = off_y;
00136         offZ = off_z;
00137       }
00138 
00139       void GeneratedMesh::parse_options(const std::vector<std::string> &groups)
00140       {
00141         for (size_t i=1; i < groups.size(); i++) {
00142           std::vector<std::string> option;
00143           stk::util::tokenize(groups[i], ":", option);
00144           // option[0] is the type of the option and option[1] is the argument to the option.
00145 
00146           if (option[0] == "shell") {
00147             // Option of the form  "shell:xXyYzZ"
00148             // The argument specifies whether there is a shell block
00149             // at the location. 'x' is minX, 'X' is maxX, etc.
00150             size_t length = option[1].size();
00151             for (size_t j=0; j < length; j++) {
00152               switch (option[1][j]) {
00153               case 'x':
00154                 add_shell_block(MX);
00155                 break;
00156               case 'X':
00157                 add_shell_block(PX);
00158                 break;
00159               case 'y':
00160                 add_shell_block(MY);
00161                 break;
00162               case 'Y':
00163                 add_shell_block(PY);
00164                 break;
00165               case 'z':
00166                 add_shell_block(MZ);
00167                 break;
00168               case 'Z':
00169                 add_shell_block(PZ);
00170                 break;
00171               default:
00172                 std::cerr << "ERROR: Unrecognized shell location option '"
00173                           << option[1][j]
00174                           << "'.";
00175               }
00176             }
00177           }
00178     else if (option[0] == "nodeset") {
00179       // Option of the form  "nodeset:xXyYzZ"
00180       // The argument specifies whether there is a nodeset
00181       // at the location. 'x' is minX, 'X' is maxX, etc.
00182       size_t length = option[1].size();
00183       for (size_t j=0; j < length; j++) {
00184         switch (option[1][j]) {
00185         case 'x':
00186     add_nodeset(MX);
00187     break;
00188         case 'X':
00189     add_nodeset(PX);
00190     break;
00191         case 'y':
00192     add_nodeset(MY);
00193     break;
00194         case 'Y':
00195     add_nodeset(PY);
00196     break;
00197         case 'z':
00198     add_nodeset(MZ);
00199     break;
00200         case 'Z':
00201     add_nodeset(PZ);
00202     break;
00203         default:
00204     std::cerr << "ERROR: Unrecognized nodeset location option '"
00205         << option[1][j]
00206         << "'.";
00207         }
00208       }
00209     }
00210     else if (option[0] == "sideset") {
00211       // Option of the form  "sideset:xXyYzZ"
00212       // The argument specifies whether there is a sideset
00213       // at the location. 'x' is minX, 'X' is maxX, etc.
00214       size_t length = option[1].size();
00215       for (size_t j=0; j < length; j++) {
00216         switch (option[1][j]) {
00217         case 'x':
00218     add_sideset(MX);
00219     break;
00220         case 'X':
00221     add_sideset(PX);
00222     break;
00223         case 'y':
00224     add_sideset(MY);
00225     break;
00226         case 'Y':
00227     add_sideset(PY);
00228     break;
00229         case 'z':
00230     add_sideset(MZ);
00231     break;
00232         case 'Z':
00233     add_sideset(PZ);
00234     break;
00235         default:
00236     std::cerr << "ERROR: Unrecognized sideset location option '"
00237         << option[1][j]
00238         << "'.";
00239         }
00240       }
00241     }
00242           else if (option[0] == "scale") {
00243             // Option of the form  "scale:xs,ys,zs
00244             std::vector<std::string> tokens;
00245             stk::util::tokenize(option[1], ",", tokens);
00246             assert(tokens.size() == 3);
00247             sclX = std::strtod(tokens[0].c_str(), NULL);
00248             sclY = std::strtod(tokens[1].c_str(), NULL);
00249             sclZ = std::strtod(tokens[2].c_str(), NULL);
00250           }
00251 
00252           else if (option[0] == "offset") {
00253             // Option of the form  "offset:xo,yo,zo
00254             std::vector<std::string> tokens;
00255             stk::util::tokenize(option[1], ",", tokens);
00256             assert(tokens.size() == 3);
00257             offX = std::strtod(tokens[0].c_str(), NULL);
00258             offY = std::strtod(tokens[1].c_str(), NULL);
00259             offZ = std::strtod(tokens[2].c_str(), NULL);
00260           }
00261 
00262           else if (option[0] == "zdecomp") {
00263             // Option of the form  "zdecomp:1,1,2,2,1,2,...
00264             // Specifies the number of intervals in the z direction
00265             // for each processor.  The number of tokens must match
00266             // the number of processors.  Note that the new numZ will
00267             // be the sum of the intervals specified in this command.
00268             std::vector<std::string> tokens;
00269             stk::util::tokenize(option[1], ",", tokens);
00270             assert(tokens.size() == processorCount);
00271             std::vector<size_t> Zs;
00272             numZ = 0;
00273             for (size_t j = 0; j < processorCount; j++) {
00274               Zs.push_back(std::strtol(tokens[j].c_str(), NULL, 10));
00275               numZ += Zs[j];
00276             }
00277             myNumZ = Zs[myProcessor];
00278             myStartZ = 0;
00279             for (size_t j=0; j < myProcessor; j++) {
00280               myStartZ += Zs[j];
00281             }
00282           }
00283 
00284           else if (option[0] == "bbox") {
00285             // Bounding-Box Option of the form  "bbox:xmin,ymin,zmin,xmax,ymax,zmaxo
00286             std::vector<std::string> tokens;
00287             stk::util::tokenize(option[1], ",", tokens);
00288             assert(tokens.size() == 6);
00289             double xmin = std::strtod(tokens[0].c_str(), NULL);
00290             double ymin = std::strtod(tokens[1].c_str(), NULL);
00291             double zmin = std::strtod(tokens[2].c_str(), NULL);
00292             double xmax = std::strtod(tokens[3].c_str(), NULL);
00293             double ymax = std::strtod(tokens[4].c_str(), NULL);
00294             double zmax = std::strtod(tokens[5].c_str(), NULL);
00295 
00296             set_bbox(xmin, ymin, zmin,  xmax, ymax, zmax);
00297           }
00298 
00299           else if (option[0] == "rotate") {
00300             // Rotate Option of the form  "rotate:axis,angle,axis,angle,...
00301             std::vector<std::string> tokens;
00302             stk::util::tokenize(option[1], ",", tokens);
00303             assert(tokens.size() %2 == 0);
00304             for (size_t ir=0; ir < tokens.size();) {
00305               std::string axis = tokens[ir++];
00306               double angle_degree = std::strtod(tokens[ir++].c_str(), NULL);
00307               set_rotation(axis, angle_degree);
00308             }
00309           }
00310 
00311           else if (option[0] == "help") {
00312       std::cerr << "\nValid Options for GeneratedMesh parameter string:\n"
00313           << "\tIxJxK -- specifies intervals; must be first option. Ex: 4x10x12\n"
00314           << "\toffset:xoff, yoff, zoff\n"
00315           << "\tscale: xscl, yscl, zscl\n"
00316           << "\tzdecomp:n1,n2,n3,...,n#proc\n"
00317           << "\tbbox: xmin, ymin, zmin, xmax, ymax, zmax\n"
00318           << "\trotate: axis,angle,axis,angle,...\n"
00319           << "\tshell:xXyYzZ (specifies which plane to apply shell)\n"
00320           << "\tnodeset:xXyXzZ (specifies which plane to apply nodeset)\n"
00321           << "\tsideset:xXyXzZ (specifies which plane to apply sideset)\n"
00322           << "\tshow -- show mesh parameters\n"
00323           << "\thelp -- show this list\n\n";
00324           }
00325 
00326           else if (option[0] == "show") {
00327       show_parameters();
00328           }
00329 
00330           else {
00331             std::cerr << "ERROR: Unrecognized option '" << option[0]
00332                       << "'.  It will be ignored.\n";
00333           }
00334         }
00335       }
00336 
00337       void GeneratedMesh::show_parameters() const
00338       {
00339   if (myProcessor == 0) {
00340     std::cerr << "\nMesh Parameters:\n"
00341         << "\tIntervals: " << numX << " by " << numY << " by " << numZ << "\n"
00342         << "\tX = " << sclX << " * (0.." << numX << ") + " << offX
00343         << "\tRange: " << offX << " <= X <= "  << offX + numX * sclX << "\n"
00344         << "\tY = " << sclY << " * (0.." << numY << ") + " << offY
00345         << "\tRange: " << offY << " <= Y <= "  << offY + numY * sclY << "\n"
00346         << "\tZ = " << sclZ << " * (0.." << numZ << ") + " << offZ
00347         << "\tRange: " << offZ << " <= Z <= "  << offZ + numZ * sclZ << "\n\n"
00348         << "\tNode Count (total)    = " << std::setw(9) << node_count() << "\n"
00349         << "\tElement Count (total) = " << std::setw(9) << element_count() << "\n"
00350         << "\tBlock Count           = " << std::setw(9) << block_count() << "\n"
00351         << "\tNodeset Count         = " << std::setw(9) << nodeset_count() << "\n"
00352         << "\tSideset Count         = " << std::setw(9) << sideset_count() << "\n\n";
00353     if (doRotation) {
00354       std::cerr << "\tRotation Matrix: \n\t" << std::scientific ;
00355       for (int ii=0; ii < 3; ii++) {
00356         for (int jj=0; jj < 3; jj++) {
00357     std::cerr << std::setw(14) << rotmat[ii][jj] << "\t";
00358         }
00359         std::cerr << "\n\t";
00360       }
00361       std::cerr << std::fixed << "\n";
00362     }
00363   }
00364       }
00365 
00366       size_t GeneratedMesh::node_count() const
00367       {
00368         return (numX+1) * (numY+1) * (numZ+1);
00369       }
00370 
00371       size_t GeneratedMesh::node_count_proc() const
00372       {
00373         return (numX+1) * (numY+1) * (myNumZ+1);
00374       }
00375 
00376       size_t GeneratedMesh::block_count() const
00377       {
00378         return shellBlocks.size() + 1;
00379       }
00380 
00381       size_t GeneratedMesh::nodeset_count() const
00382       {
00383   return nodesets.size();
00384       }
00385       
00386       size_t GeneratedMesh::sideset_count() const
00387       {
00388   return sidesets.size();
00389       }
00390 
00391       size_t GeneratedMesh::element_count() const
00392       {
00393         size_t count = element_count(1);
00394         for (size_t i=0; i < shellBlocks.size(); i++) {
00395           count += element_count(i+2);
00396         }
00397         return count;
00398       }
00399 
00400       size_t GeneratedMesh::element_count_proc() const
00401       {
00402         size_t count = 0;
00403         for (size_t i=0; i < block_count(); i++) {
00404           count += element_count_proc(i+1);
00405         }
00406         return count;
00407       }
00408 
00409       size_t GeneratedMesh::element_count(size_t block_number) const
00410       {
00411         assert(block_number <= block_count());
00412 
00413         if (block_number == 1) {
00414           return numX * numY * numZ;
00415         } else {
00416           ShellLocation loc = shellBlocks[block_number-2];
00417           return shell_element_count(loc);
00418         }
00419       }
00420 
00421       size_t GeneratedMesh::shell_element_count(ShellLocation loc) const
00422       {
00423         switch (loc) {
00424   case MX:
00425   case PX:
00426     return numY * numZ;
00427   case MY:
00428   case PY:
00429     return numX * numZ;
00430   case MZ:
00431   case PZ:
00432     return numX * numY;
00433         }
00434         return 0;
00435       }
00436 
00437       size_t GeneratedMesh::element_count_proc(size_t block_number) const
00438       {
00439         assert(block_number <= block_count());
00440 
00441         if (block_number == 1) {
00442           return numX * numY * myNumZ;
00443         } else {
00444           ShellLocation loc = shellBlocks[block_number-2];
00445           return shell_element_count_proc(loc);
00446         }
00447       }
00448 
00449       size_t GeneratedMesh::shell_element_count_proc(ShellLocation loc) const
00450       {
00451         switch (loc) {
00452   case MX:
00453   case PX:
00454     return numY * myNumZ;
00455   case MY:
00456   case PY:
00457     return numX * myNumZ;
00458   case MZ:
00459     if (myProcessor == 0)
00460       return numX * numY;
00461     else
00462       return 0;
00463   case PZ:
00464     if (myProcessor == processorCount -1)
00465       return numX * numY;
00466     else
00467       return 0;
00468         }
00469         return 0;
00470       }
00471 
00472       size_t GeneratedMesh::nodeset_node_count(size_t id) const
00473       {
00474   // id is position in nodeset list + 1
00475   assert(id > 0 && id <= nodesets.size());
00476   ShellLocation loc = nodesets[id-1];
00477   switch (loc) {
00478   case MX:
00479   case PX:
00480     return (numY+1) * (numZ+1);
00481   case MY:
00482   case PY:
00483     return (numX+1) * (numZ+1);
00484   case MZ:
00485   case PZ:
00486     return (numX+1) * (numY+1);
00487   }
00488   return 0;
00489       }
00490 
00491       size_t GeneratedMesh::nodeset_node_count_proc(size_t id) const
00492       {
00493   // id is position in nodeset list + 1
00494   assert(id > 0 && id <= nodesets.size());
00495   ShellLocation loc = nodesets[id-1];
00496   switch (loc) {
00497   case MX:
00498   case PX:
00499     return (numY+1) * (myNumZ+1);
00500   case MY:
00501   case PY:
00502     return (numX+1) * (myNumZ+1);
00503   case MZ:
00504     if (myProcessor == 0)
00505       return (numX+1) * (numY+1);
00506     else
00507       return 0;
00508   case PZ:
00509     if (myProcessor == processorCount -1)
00510       return (numX+1) * (numY+1);
00511     else
00512       return 0;
00513   }
00514   return 0;
00515       }
00516 
00517       size_t GeneratedMesh::sideset_side_count(size_t id) const
00518       {
00519   // id is position in sideset list + 1
00520   assert(id > 0 && id <= sidesets.size());
00521   ShellLocation loc = sidesets[id-1];
00522   switch (loc) {
00523   case MX:
00524   case PX:
00525     return numY * numZ;
00526   case MY:
00527   case PY:
00528     return numX * numZ;
00529   case MZ:
00530   case PZ:
00531     return numX * numY;
00532   }
00533   return 0;
00534       }
00535 
00536       size_t GeneratedMesh::sideset_side_count_proc(size_t id) const
00537       {
00538   // id is position in sideset list + 1
00539   assert(id > 0 && id <= sidesets.size());
00540   ShellLocation loc = sidesets[id-1];
00541   switch (loc) {
00542   case MX:
00543   case PX:
00544     return numY * myNumZ;
00545   case MY:
00546   case PY:
00547     return numX * myNumZ;
00548   case MZ:
00549     if (myProcessor == 0)
00550       return numX * numY;
00551     else
00552       return 0;
00553   case PZ:
00554     if (myProcessor == processorCount -1)
00555       return numX * numY;
00556     else
00557       return 0;
00558   }
00559   return 0;
00560       }
00561 
00562       std::pair<std::string,int> GeneratedMesh::topology_type(size_t block_number) const
00563       {
00564         assert(block_number <= block_count() && block_number > 0);
00565 
00566         if (block_number == 1) {
00567           return std::make_pair(std::string("hex8"), 8);;
00568         } else {
00569           return std::make_pair(std::string("shell4"), 4);
00570         }
00571       }
00572 
00573       void GeneratedMesh::node_map(std::vector<int> &map)
00574       {
00575         size_t count = node_count_proc();
00576         map.resize(count);
00577         size_t offset = myStartZ * (numX+1) * (numY+1);
00578         for (size_t i=0; i < count; i++) {
00579           map[i] = static_cast<int>(offset + i + 1);
00580         }
00581       }
00582 
00583       size_t GeneratedMesh::communication_node_count_proc() const
00584       {
00585   size_t count = (numX+1) * (numY+1);
00586   if (myProcessor != 0 && myProcessor != processorCount-1)
00587     count *= 2;
00588     
00589   return count;
00590       }
00591 
00592       void GeneratedMesh::node_communication_map(std::vector<int> &map, std::vector<int> &proc)
00593       {
00594         size_t count = (numX+1) * (numY+1);
00595         size_t slab = count;
00596         if (myProcessor != 0 && myProcessor != processorCount-1)
00597           count *= 2;
00598 
00599         map.resize(count);
00600         proc.resize(count);
00601         size_t j = 0;
00602         if (myProcessor != 0) {
00603           size_t offset = myStartZ * (numX+1) * (numY+1);
00604           for (size_t i=0; i < slab; i++) {
00605             map[j] = static_cast<int>(offset + i + 1);
00606             proc[j++] = static_cast<int>(myProcessor-1);
00607           }
00608         }
00609         if (myProcessor != processorCount-1) {
00610           size_t offset = (myStartZ + myNumZ) * (numX+1) * (numY+1);
00611           for (size_t i=0; i < slab; i++) {
00612             map[j] = static_cast<int>(offset + i + 1);
00613             proc[j++] = static_cast<int>(myProcessor+1);
00614           }
00615         }
00616       }
00617 
00618       void GeneratedMesh::element_map(size_t block_number, std::vector<int> &map) const
00619       {
00620         assert(block_number <= block_count() && block_number > 0);
00621 
00622         size_t count = element_count_proc(block_number);
00623         map.resize(count);
00624 
00625         if (block_number == 1) {
00626           // Hex block...
00627           count = element_count_proc(1);
00628           size_t offset = myStartZ * numX * numY;
00629           for (size_t i=0; i < count; i++) {
00630             map[i] = static_cast<int>(offset + i + 1);
00631           }
00632         } else {
00633           size_t start = element_count(1);
00634 
00635           // Shell blocks...
00636           for (size_t ib=0; ib < shellBlocks.size(); ib++) {
00637             count = element_count_proc(ib+2);
00638             if (static_cast<size_t>(block_number) == ib + 2) {
00639               size_t offset = 0;
00640               ShellLocation loc = shellBlocks[ib];
00641 
00642               switch (loc) {
00643               case MX:
00644               case PX:
00645                 offset = myStartZ * numY;
00646                 break;
00647 
00648               case MY:
00649               case PY:
00650                 offset = myStartZ * numX;
00651                 break;
00652 
00653               case MZ:
00654               case PZ:
00655                 offset = 0;
00656                 break;
00657               }
00658               for (size_t i=0; i < count; i++) {
00659                 map[i] = static_cast<int>(start + offset + i + 1);
00660               }
00661             } else {
00662               start += element_count(ib+2);
00663             }
00664           }
00665         }
00666       }
00667 
00668       void GeneratedMesh::element_map(std::vector<int> &map) const
00669       {
00670         size_t count = element_count_proc();
00671         map.resize(count);
00672 
00673         size_t k = 0;
00674         // Hex block...
00675         count = element_count_proc(1);
00676         size_t offset = myStartZ * numX * numY;
00677         for (size_t i=0; i < count; i++) {
00678           map[k++] = static_cast<int>(offset + i + 1);
00679         }
00680 
00681         size_t start = element_count(1);
00682 
00683         // Shell blocks...
00684         for (size_t ib=0; ib < shellBlocks.size(); ib++) {
00685           count = element_count_proc(ib+2);
00686           offset = 0;
00687           ShellLocation loc = shellBlocks[ib];
00688 
00689           switch (loc) {
00690           case MX:
00691           case PX:
00692             offset = myStartZ * numY;
00693             break;
00694 
00695           case MY:
00696           case PY:
00697             offset = myStartZ * numX;
00698             break;
00699 
00700           case MZ:
00701           case PZ:
00702             offset = 0;
00703             break;
00704           }
00705           for (size_t i=0; i < count; i++) {
00706             map[k++] = static_cast<int>(start + offset + i + 1);
00707           }
00708           start += element_count(ib+2);
00709         }
00710       }
00711 
00712       void GeneratedMesh::element_surface_map(ShellLocation loc, std::vector<int> &map) const
00713       {
00714   size_t count = shell_element_count_proc(loc);
00715   map.resize(2*count);
00716 
00717   size_t index  = 0;
00718   size_t offset = 0;
00719 
00720   switch (loc) {
00721   case MX:
00722     offset = myStartZ * numX * numY + 1;  // 1-based elem id
00723     for (size_t k = 0; k < myNumZ; ++k) {
00724       for (size_t j = 0; j < numY; ++j) {
00725         map[index++] = offset; 
00726         map[index++] = 3; // 0-based local face id
00727         offset += numX;
00728       }
00729     }
00730     break;
00731 
00732   case PX:
00733     offset = myStartZ * numX * numY + numX;
00734     for (size_t k = 0; k < myNumZ; ++k) {
00735       for (size_t j = 0; j < numY; ++j) {
00736         map[index++] = offset; // 1-based elem id
00737         map[index++] = 1; // 0-based local face id
00738         offset += numX;
00739       }
00740     }
00741     break;
00742       
00743   case MY:
00744     offset = myStartZ * numX * numY + 1;
00745     for (size_t k = 0; k < myNumZ; ++k) {
00746       for (size_t i = 0; i < numX; ++i) {
00747         map[index++] = offset++;
00748         map[index++] = 0; // 0-based local face id
00749       }
00750       offset+= numX * (numY-1);
00751     }
00752     break;
00753 
00754   case PY:
00755     offset = myStartZ * numX * numY + numX * (numY-1) +1;
00756     for (size_t k = 0; k < myNumZ; ++k) {
00757       for (size_t i = 0; i < numX; ++i) {
00758         map[index++] = offset++;
00759         map[index++] = 2; // 0-based local face id
00760       }
00761       offset+= numX * (numY-1);
00762     }
00763     break;
00764 
00765   case MZ:
00766     if (myProcessor == 0) {
00767       offset = 1;
00768       for (size_t i=0; i < numY; i++) {
00769         for (size_t j=0; j < numX; j++) {
00770     map[index++] = offset++;
00771     map[index++] = 4;
00772         }
00773       }
00774     }
00775     break;
00776 
00777   case PZ:
00778     if (myProcessor == processorCount-1) {
00779       offset = (numZ-1)*numX*numY + 1;
00780       for (size_t i=0, k=0; i < numY; i++) {
00781         for (size_t j=0; j < numX; j++, k++) {
00782     map[index++] = offset++;
00783     map[index++] = 5;
00784         }
00785       }
00786     }
00787     break;
00788   }
00789       }
00790 
00791       void GeneratedMesh::coordinates(std::vector<double> &coord) const
00792       {
00793         /* create global coordinates */
00794         size_t count = node_count_proc();
00795         coord.resize(count * 3);
00796 
00797         size_t k = 0;
00798         for (size_t m=myStartZ; m < myStartZ+myNumZ+1; m++) {
00799           for (size_t i=0; i < numY+1; i++) {
00800             for (size_t j=0; j < numX+1; j++) {
00801               coord[k++] = sclX * static_cast<double>(j) + offX;
00802               coord[k++] = sclY * static_cast<double>(i) + offY;
00803               coord[k++] = sclZ * static_cast<double>(m) + offZ;
00804             }
00805           }
00806         }
00807 
00808         if (doRotation) {
00809           for (size_t i=0; i < count*3; i+=3) {
00810             double xn = coord[i+0];
00811             double yn = coord[i+1];
00812             double zn = coord[i+2];
00813             coord[i+0] = xn * rotmat[0][0] + yn * rotmat[1][0] + zn * rotmat[2][0];
00814             coord[i+1] = xn * rotmat[0][1] + yn * rotmat[1][1] + zn * rotmat[2][1];
00815             coord[i+2] = xn * rotmat[0][2] + yn * rotmat[1][2] + zn * rotmat[2][2];
00816           }
00817         }
00818       }
00819 
00820       void GeneratedMesh::coordinates(std::vector<double> &x,
00821                                       std::vector<double> &y,
00822                                       std::vector<double> &z) const
00823       {
00824         /* create global coordinates */
00825         size_t count = node_count_proc();
00826         x.resize(count);
00827         y.resize(count);
00828         z.resize(count);
00829 
00830         size_t k = 0;
00831         for (size_t m=myStartZ; m < myStartZ+myNumZ+1; m++) {
00832           for (size_t i=0; i < numY+1; i++) {
00833             for (size_t j=0; j < numX+1; j++) {
00834               x[k] = sclX * static_cast<double>(j) + offX;
00835               y[k] = sclY * static_cast<double>(i) + offY;
00836               z[k] = sclZ * static_cast<double>(m) + offZ;
00837               ++k;
00838             }
00839           }
00840         }
00841         if (doRotation) {
00842           for (size_t i=0; i < count; i++) {
00843             double xn = x[i];
00844             double yn = y[i];
00845             double zn = z[i];
00846             x[i] = xn * rotmat[0][0] + yn * rotmat[1][0] + zn * rotmat[2][0];
00847             y[i] = xn * rotmat[0][1] + yn * rotmat[1][1] + zn * rotmat[2][1];
00848             z[i] = xn * rotmat[0][2] + yn * rotmat[1][2] + zn * rotmat[2][2];
00849           }
00850         }
00851       }
00852 
00853       void GeneratedMesh::connectivity(size_t block_number, std::vector<int> &connect) const
00854       {
00855   assert(block_number <= block_count());
00856 
00857   size_t xp1yp1 = (numX+1) * (numY+1);
00858 
00859   /* build connectivity array (node list) for mesh */
00860   if (block_number == 1) {  // HEX Element Block
00861     connect.resize(element_count_proc(block_number)*8);
00862 
00863     size_t cnt = 0;
00864     for (size_t m=myStartZ; m < myNumZ+myStartZ; m++) {
00865       for (size_t i=0, k=0; i < numY; i++) {
00866         for (size_t j=0; j < numX; j++, k++) {
00867     size_t base = (m*xp1yp1) + k + i + 1;
00868     ;
00869     connect[cnt++] = static_cast<int>(base);
00870     connect[cnt++] = static_cast<int>(base+1);
00871     connect[cnt++] = static_cast<int>(base+numX+2);
00872     connect[cnt++] = static_cast<int>(base+numX+1);
00873 
00874     connect[cnt++] = static_cast<int>(xp1yp1 + base);
00875     connect[cnt++] = static_cast<int>(xp1yp1 + base+1);
00876     connect[cnt++] = static_cast<int>(xp1yp1 + base+numX+2);
00877     connect[cnt++] = static_cast<int>(xp1yp1 + base+numX+1);
00878         }
00879       }
00880     }
00881   } else { // Shell blocks....
00882     ShellLocation loc = shellBlocks[block_number-2];
00883     connect.resize(element_count_proc(block_number)*4);
00884 
00885     size_t cnt = 0;
00886     switch (loc) {
00887     case MX:  // Minumum X Face
00888       for (size_t i=0; i < myNumZ; i++) {
00889         size_t layer_off = i * xp1yp1;
00890         for (size_t j=0; j < numY; j++) {
00891     size_t base = layer_off + j * (numX+1) + 1 + myStartZ * xp1yp1;
00892     connect[cnt++] = static_cast<int>(base);
00893     connect[cnt++] = static_cast<int>(base + xp1yp1);
00894     connect[cnt++] = static_cast<int>(base + xp1yp1 + (numX+1));
00895     connect[cnt++] = static_cast<int>(base + (numX+1));
00896         }
00897       }
00898       break;
00899     case PX: // Maximum X Face
00900       for (size_t i=0; i < myNumZ; i++) {
00901         size_t layer_off = i * xp1yp1;
00902         for (size_t j=0; j < numY; j++) {
00903     size_t base = layer_off + j * (numX+1) + numX + 1 + myStartZ * xp1yp1;
00904     connect[cnt++] = static_cast<int>(base);
00905     connect[cnt++] = static_cast<int>(base + (numX+1));
00906     connect[cnt++] = static_cast<int>(base + xp1yp1 + (numX+1));
00907     connect[cnt++] = static_cast<int>(base + xp1yp1);
00908         }
00909       }
00910       break;
00911     case MY: // Minumum Y Face
00912       for (size_t i=0; i < myNumZ; i++) {
00913         size_t layer_off = i * xp1yp1;
00914         for (size_t j=0; j < numX; j++) {
00915     size_t base = layer_off + j + 1 + myStartZ * xp1yp1;
00916     connect[cnt++] = static_cast<int>(base);
00917     connect[cnt++] = static_cast<int>(base + 1);
00918     connect[cnt++] = static_cast<int>(base + xp1yp1 + 1);
00919     connect[cnt++] = static_cast<int>(base + xp1yp1);
00920         }
00921       }
00922       break;
00923     case PY: // Maximum Y Face
00924       for (size_t i=0; i < myNumZ; i++) {
00925         size_t layer_off = i * xp1yp1;
00926         for (size_t j=0; j < numX; j++) {
00927     size_t base = layer_off + (numX+1)*(numY) + j + 1 + myStartZ * xp1yp1;
00928     connect[cnt++] = static_cast<int>(base);
00929     connect[cnt++] = static_cast<int>(base + xp1yp1);
00930     connect[cnt++] = static_cast<int>(base + xp1yp1 + 1);
00931     connect[cnt++] = static_cast<int>(base + 1);
00932         }
00933       }
00934       break;
00935     case MZ: // Minumum Z Face
00936       if (myProcessor == 0) {
00937         for (size_t i=0, k=0; i < numY; i++) {
00938     for (size_t j=0; j < numX; j++, k++) {
00939       size_t base = i + k + 1 + myStartZ * xp1yp1;
00940       connect[cnt++] = static_cast<int>(base);
00941       connect[cnt++] = static_cast<int>(base+numX+1);
00942       connect[cnt++] = static_cast<int>(base+numX+2);
00943       connect[cnt++] = static_cast<int>(base+1);
00944     }
00945         }
00946       }
00947       break;
00948     case PZ: // Maximum Z Face
00949       if (myProcessor == processorCount-1) {
00950         for (size_t i=0, k=0; i < numY; i++) {
00951     for (size_t j=0; j < numX; j++, k++) {
00952       size_t base = xp1yp1 * (numZ - myStartZ) + k + i + 1 + myStartZ * xp1yp1;
00953       connect[cnt++] = static_cast<int>(base);
00954       connect[cnt++] = static_cast<int>(base+1);
00955       connect[cnt++] = static_cast<int>(base+numX+2);
00956       connect[cnt++] = static_cast<int>(base+numX+1);
00957     }
00958         }
00959       }
00960       break;
00961     }
00962     assert(cnt == 4 * element_count_proc(block_number));
00963   }
00964   return;
00965       }
00966 
00967       void GeneratedMesh::nodeset_nodes(size_t id, std::vector<int> &nodes) const
00968       {
00969   // id is position in nodeset list + 1
00970   assert(id > 0 && id <= nodesets.size());
00971   ShellLocation loc = nodesets[id-1];
00972   nodes.resize(nodeset_node_count_proc(id));
00973 
00974   size_t xp1yp1 = (numX+1) * (numY+1);
00975   size_t k = 0;
00976 
00977   switch (loc) {
00978   case MX:  // Minumum X Face
00979     for (size_t i=0; i < myNumZ+1; i++) {
00980       size_t layer_off = myStartZ * xp1yp1 + i * xp1yp1;
00981       for (size_t j=0; j < numY+1; j++) {
00982         nodes[k++] = layer_off + j * (numX+1) + 1;
00983       }
00984     }
00985     break;
00986   case PX: // Maximum X Face
00987     for (size_t i=0; i < myNumZ+1; i++) {
00988       size_t layer_off = myStartZ * xp1yp1 + i * xp1yp1;
00989       for (size_t j=0; j < numY+1; j++) {
00990         nodes[k++] = layer_off + j * (numX+1) + numX + 1;
00991       }
00992     }
00993     break;
00994   case MY: // Minumum Y Face
00995     for (size_t i=0; i < myNumZ+1; i++) {
00996       size_t layer_off = myStartZ * xp1yp1 + i * xp1yp1;
00997       for (size_t j=0; j < numX+1; j++) {
00998         nodes[k++] = layer_off + j + 1;
00999       }
01000     }
01001     break;
01002   case PY: // Maximum Y Face
01003     for (size_t i=0; i < myNumZ+1; i++) {
01004       size_t layer_off = myStartZ * xp1yp1 + i * xp1yp1;
01005       for (size_t j=0; j < numX+1; j++) {
01006         nodes[k++] = layer_off + (numX+1)*(numY) + j + 1;
01007       }
01008     }
01009     break;
01010   case MZ: // Minumum Z Face
01011     if (myProcessor == 0) {
01012       for (size_t i=0; i < (numY+1) * (numX+1); i++) {
01013         nodes[i] = i+1;
01014       }
01015     }
01016     break;
01017   case PZ: // Maximum Z Face
01018     if (myProcessor == processorCount-1) {
01019       size_t offset = (numY+1) * (numX+1) * numZ;
01020       for (size_t i=0; i < (numY+1) * (numX+1); i++) {
01021         nodes[i] = offset + i+1;
01022       }
01023     }
01024     break;
01025   }
01026       }
01027 
01028       void GeneratedMesh::sideset_elem_sides(size_t id, std::vector<int> &elem_sides) const
01029       {
01030   // id is position in sideset list + 1
01031   assert(id > 0 && id <= sidesets.size());
01032   ShellLocation loc = sidesets[id-1];
01033     
01034   // If there is a shell block on this face, then the sideset is
01035   // applied to the shell block; if not, it is applied to the
01036   // underlying hex elements.
01037   bool underlying_shell = false;
01038   size_t shell_block = 0;
01039   for (size_t i = 0; i < shellBlocks.size(); i++) {
01040     if (shellBlocks[i] == loc) {
01041       underlying_shell = true;
01042       shell_block = i+2;
01043       break;
01044     }
01045   }
01046 
01047   if (underlying_shell) {
01048     // Get ids of shell elements at this location...
01049     element_map(shell_block, elem_sides);
01050       
01051     // Insert face_ordinal in between each entry in elem_sides...
01052     // Face will be 0 for all shells...
01053     elem_sides.resize(2*sideset_side_count_proc(id));
01054     int face_ordinal = 0;
01055     int i = 2* (int)sideset_side_count_proc(id) - 1;
01056     int j =    (int)sideset_side_count_proc(id) - 1;
01057     while (i >= 0) {
01058       elem_sides[i--] = face_ordinal;
01059       elem_sides[i--] = elem_sides[j--];
01060     }
01061   } else {
01062     element_surface_map(loc, elem_sides);
01063   }
01064       }
01065 
01066       void GeneratedMesh::set_rotation(const std::string &axis, double angle_degrees)
01067       {
01068         // PI / 180. Used in converting angle in degrees to radians
01069         static double degang = std::atan2(0.0, -1.0) / 180.0;
01070 
01071         doRotation = true;
01072 
01073         int n1 = -1;
01074   int n2 = -1;
01075   int n3 = -1;
01076 
01077         if (axis == "x" || axis == "X") {
01078           n1 = 1; n2 = 2; n3 = 0;
01079         } else if (axis == "y" || axis == "Y") {
01080           n1 = 2; n2 = 0; n3 = 1;
01081         } else if (axis == "z" || axis == "Z") {
01082           n1 = 0; n2 = 1; n3 = 2;
01083         } else {
01084           std::cerr << "\nInvalid axis specification '" << axis << "'. Valid options are 'x', 'y', or 'z'\n";
01085           return;
01086         }
01087 
01088         double ang = angle_degrees * degang; // Convert angle in degrees to radians
01089         double cosang = std::cos(ang);
01090         double sinang = std::sin(ang);
01091 
01092   assert(n1 >= 0 && n2 >= 0 && n3 >= 0);
01093         double by[3][3];
01094         by[n1][n1] =  cosang;
01095         by[n2][n1] = -sinang;
01096         by[n1][n3] = 0.0;
01097         by[n1][n2] =  sinang;
01098         by[n2][n2] =  cosang;
01099         by[n2][n3] = 0.0;
01100         by[n3][n1] = 0.0;
01101         by[n3][n2] = 0.0;
01102         by[n3][n3] = 1.0;
01103 
01104         double res[3][3];
01105         for (size_t i=0; i < 3; i++) {
01106           res[i][0] = rotmat[i][0]*by[0][0] + rotmat[i][1]*by[1][0] + rotmat[i][2]*by[2][0];
01107           res[i][1] = rotmat[i][0]*by[0][1] + rotmat[i][1]*by[1][1] + rotmat[i][2]*by[2][1];
01108           res[i][2] = rotmat[i][0]*by[0][2] + rotmat[i][1]*by[1][2] + rotmat[i][2]*by[2][2];
01109         }
01110 
01111 #if 1
01112         std::memcpy(rotmat, res, 9*sizeof(double));
01113 #else
01114         for (int i=0; i < 3; i++) {
01115           for (int j=0; j < 3; j++) {
01116             rotmat[i][j] = res[i][j];
01117           }
01118         }
01119 #endif
01120       }
01121     }
01122   }
01123 }
01124 
01125 #if defined(DEBUG)
01126 #include <sstream>
01127   std::string ToString(size_t t) {
01128     std::ostringstream os;
01129     os << t;
01130     return os.str();
01131   }
01132 
01133 #include <exodusII.h>
01134 #include <ne_nemesisI.h>
01135 
01136   int main() {
01137     int num_processors = 8;
01138     for (int proc = 0; proc < num_processors; proc++) {
01139 
01140       stk::GeneratedMesh mesh(100, 125, 10*num_processors, num_processors, proc);
01141 
01142       std::cerr << "Node Count (total)    = " << mesh.node_count() << "\n";
01143       std::cerr << "Node Count (proc)     = " << mesh.node_count_proc() << "\n";
01144       std::cerr << "Element Count (total) = " << mesh.element_count() << "\n";
01145       std::cerr << "Element Count (proc)  = " << mesh.element_count_proc() << "\n";
01146       std::cerr << "Block Count           = " << mesh.block_count() << "\n";
01147 
01148       int CPU_word_size = 8;                   /* sizeof(float) */
01149       int IO_word_size = 8;                    /* (4 bytes) */
01150       std::string name = "test-scale.e";
01151       if (num_processors > 1) {
01152         name += "." + ToString(num_processors) + "." + ToString(proc);
01153       }
01154       int exoid = ex_create (name.c_str(), EX_CLOBBER, &CPU_word_size, &IO_word_size);
01155 
01156       int num_nodes = mesh.node_count_proc();
01157       int num_elems = mesh.element_count_proc();
01158       int num_elem_blk = mesh.block_count();
01159       int error = ex_put_init (exoid, "title", 3, num_nodes, num_elems,
01160                                num_elem_blk, 0, 0);
01161 
01162       if (num_processors > 1) {
01163         std::vector<int> nodes;
01164         std::vector<int> procs;
01165         mesh.node_communication_map(nodes, procs);
01166 
01167         int node_map_ids[1] = {1};
01168         int node_map_node_cnts[1] = {procs.size()};
01169         ne_put_init_info(exoid, num_processors, 1, "p");
01170         ne_put_loadbal_param(exoid, 0, 0, 0, 0, 0, 1, 0, proc);
01171         ne_put_cmap_params(exoid, node_map_ids, node_map_node_cnts, 0, 0, proc);
01172         ne_put_node_cmap(exoid, 1, &nodes[0], &procs[0], proc);
01173       }
01174 
01175       for (int i=1; i < mesh.block_count(); i++) {
01176         std::cerr << "Block " << i+1 << " has " << mesh.element_count_proc(i+1) << " "
01177                   << mesh.topology_type(i+1).first <<" elements\n";
01178 
01179         std::string btype = mesh.topology_type(i+1).first;
01180         error = ex_put_elem_block(exoid, i+1, btype.c_str(),
01181                                   mesh.element_count_proc(i+1),
01182                                   mesh.topology_type(i+1).second, 0);
01183       }
01184       {
01185         std::cerr << "Block " << 1 << " has " << mesh.element_count_proc(1) << " "
01186                   << mesh.topology_type(1).first <<" elements\n";
01187 
01188         std::string btype = mesh.topology_type(1).first;
01189         error = ex_put_elem_block(exoid, 1, btype.c_str(),
01190                                   mesh.element_count_proc(1),
01191                                   mesh.topology_type(1).second, 0);
01192       }
01193 
01194       if (num_processors > 1) {
01195         std::vector<int> map;
01196         mesh.node_map(map);
01197         ex_put_id_map(exoid, EX_NODE_MAP, &map[0]);
01198 
01199         mesh.element_map(map);
01200         ex_put_id_map(exoid, EX_ELEM_MAP, &map[0]);
01201       }
01202 
01203       std::cerr << "Outputting connectivity...\n";
01204       for (int i=1; i < mesh.block_count(); i++) {
01205         if (mesh.element_count_proc(i+1) > 0) {
01206           std::vector<int> connectivity;
01207           mesh.connectivity(i+1, connectivity);
01208           ex_put_elem_conn(exoid, i+1, &connectivity[0]);
01209         }
01210       }
01211       {
01212         std::vector<int> connectivity;
01213         mesh.connectivity(1, connectivity);
01214         ex_put_elem_conn(exoid, 1, &connectivity[0]);
01215       }
01216 
01217       {
01218         std::vector<double> x;
01219         std::vector<double> y;
01220         std::vector<double> z;
01221         mesh.coordinates(x, y, z);
01222         error = ex_put_coord (exoid, &x[0], &y[0], &z[0]);
01223       }
01224 
01225       ex_close(exoid);
01226     }
01227   }
01228 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends