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/unit_test_support/tokenize.hpp>
00011 
00012 #include <cmath>
00013 #include <cstring>
00014 
00015 #include <vector>
00016 #include <string>
00017 #include <iostream>
00018 #include <iomanip>
00019 #include <assert.h>
00020 
00021 namespace stk {
00022   namespace io {
00023     namespace util {
00024       GeneratedMesh::GeneratedMesh(int num_x, int num_y, int num_z,
00025                                    int proc_count, int my_proc) :
00026         numX(num_x), numY(num_y), numZ(num_z), myNumZ(num_z), myStartZ(0),
00027         processorCount(proc_count), myProcessor(my_proc),
00028         offX(0), offY(0), offZ(0),
00029         sclX(1), sclY(1), sclZ(1),
00030         doRotation(false)
00031       {
00032         initialize();
00033       }
00034 
00035       GeneratedMesh::GeneratedMesh(const std::string &parameters,
00036                                    int proc_count, int my_proc) :
00037         numX(0), numY(0), numZ(0), myNumZ(0), myStartZ(0),
00038         processorCount(proc_count), myProcessor(my_proc),
00039         offX(0), offY(0), offZ(0),
00040         sclX(1), sclY(1), sclZ(1),
00041         doRotation(false)
00042       {
00043         std::vector<std::string> groups = stk::io::util::tokenize(parameters, stk::io::util::recognize("|+"));
00044 
00045         // First 'group' is the interval specification -- IxJxK
00046         std::vector<std::string> tokens = stk::io::util::tokenize(groups[0], stk::io::util::recognize("x"));
00047         assert(tokens.size() == 3);
00048         numX = std::strtol(tokens[0].c_str(), NULL, 10);
00049         numY = std::strtol(tokens[1].c_str(), NULL, 10);
00050         numZ = std::strtol(tokens[2].c_str(), NULL, 10);
00051 
00052         initialize();
00053         parse_options(groups);
00054 
00055       }
00056 
00057       GeneratedMesh::~GeneratedMesh() {}
00058 
00059       void GeneratedMesh::initialize()
00060       {
00061         assert(numZ >= processorCount);
00062         if (processorCount > 1) {
00063           myNumZ = numZ / processorCount;
00064           if (myProcessor < (numZ % processorCount)) myNumZ++;
00065 
00066           // Determine myStartZ for this processor...
00067           size_t extra = numZ % processorCount;
00068           if (extra > myProcessor)
00069             extra = myProcessor;
00070           size_t per_proc  = numZ / processorCount;
00071           myStartZ = myProcessor * per_proc + extra;
00072         } else {
00073           myNumZ = numZ;
00074         }
00075 
00076         for (int i=0; i < 3; i++) {
00077           for (int j=0; j < 3; j++) {
00078             rotmat[i][j] = 0.0;
00079           }
00080           rotmat[i][i] = 1.0;
00081         }
00082 
00083         faceNodes[0][0] = 1;
00084         faceNodes[0][1] = 2;
00085         faceNodes[0][2] = 6;
00086         faceNodes[0][3] = 5;
00087 
00088         faceNodes[1][0] = 2;
00089         faceNodes[1][1] = 3;
00090         faceNodes[1][2] = 7;
00091         faceNodes[1][3] = 6;
00092 
00093         faceNodes[2][0] = 3;
00094         faceNodes[2][1] = 4;
00095         faceNodes[2][2] = 8;
00096         faceNodes[2][3] = 7;
00097 
00098         faceNodes[3][0] = 4;
00099         faceNodes[3][1] = 1;
00100         faceNodes[3][2] = 5;
00101         faceNodes[3][3] = 8;
00102 
00103         faceNodes[4][0] = 4;
00104         faceNodes[4][1] = 3;
00105         faceNodes[4][2] = 2;
00106         faceNodes[4][3] = 1;
00107 
00108         faceNodes[5][0] = 5;
00109         faceNodes[5][1] = 6;
00110         faceNodes[5][2] = 7;
00111         faceNodes[5][3] = 8;
00112       }
00113 
00114       size_t GeneratedMesh::add_shell_block(ShellLocation loc)
00115       {
00116         shellBlocks.push_back(loc);
00117         return shellBlocks.size();
00118       }
00119 
00120       void GeneratedMesh::set_bbox(double xmin, double ymin, double zmin,
00121                                    double xmax, double ymax, double zmax)
00122       {
00123         // NOTE: All calculations are based on the currently
00124         // active interval settings. If scale or offset or zdecomp
00125         // specified later in the option list, you may not get the
00126         // desired bounding box.
00127         double x_range = xmax - xmin;
00128         double y_range = ymax - ymin;
00129         double z_range = zmax - zmin;
00130 
00131         sclX = x_range / static_cast<double>(numX);
00132         sclY = y_range / static_cast<double>(numY);
00133         sclZ = z_range / static_cast<double>(numZ);
00134 
00135         offX = xmin;
00136         offY = ymin;
00137         offZ = zmin;
00138       }
00139 
00140       void GeneratedMesh::set_scale(double scl_x, double scl_y, double scl_z)
00141       {
00142         sclX = scl_x;
00143         sclY = scl_y;
00144         sclZ = scl_z;
00145       }
00146 
00147       void GeneratedMesh::set_offset(double off_x, double off_y, double off_z)
00148       {
00149         offX = off_x;
00150         offY = off_y;
00151         offZ = off_z;
00152       }
00153 
00154       void GeneratedMesh::parse_options(const std::vector<std::string> &groups)
00155       {
00156         for (size_t i=1; i < groups.size(); i++) {
00157           std::vector<std::string> option = stk::io::util::tokenize(groups[i], stk::io::util::recognize(":"));
00158           // option[0] is the type of the option and option[1] is the argument to the option.
00159 
00160           if (option[0] == "shell") {
00161             // Option of the form  "shell:xXyYzZ"
00162             // The argument specifies whether there is a shell block
00163             // at the location. 'x' is minX, 'X' is maxX, etc.
00164             size_t length = option[1].size();
00165             for (size_t j=0; j < length; j++) {
00166               switch (option[1][j]) {
00167               case 'x':
00168                 add_shell_block(MX);
00169                 break;
00170               case 'X':
00171                 add_shell_block(PX);
00172                 break;
00173               case 'y':
00174                 add_shell_block(MY);
00175                 break;
00176               case 'Y':
00177                 add_shell_block(PY);
00178                 break;
00179               case 'z':
00180                 add_shell_block(MZ);
00181                 break;
00182               case 'Z':
00183                 add_shell_block(PZ);
00184                 break;
00185               default:
00186                 std::cerr << "ERROR: Unrecognized shell location option '"
00187                           << option[1][j]
00188                           << "'.";
00189               }
00190             }
00191           }
00192           else if (option[0] == "scale") {
00193             // Option of the form  "scale:xs,ys,zs
00194             std::vector<std::string> tokens = stk::io::util::tokenize(option[1], stk::io::util::recognize(","));
00195             assert(tokens.size() == 3);
00196             sclX = std::strtod(tokens[0].c_str(), NULL);
00197             sclY = std::strtod(tokens[1].c_str(), NULL);
00198             sclZ = std::strtod(tokens[2].c_str(), NULL);
00199           }
00200 
00201           else if (option[0] == "offset") {
00202             // Option of the form  "offset:xo,yo,zo
00203             std::vector<std::string> tokens = stk::io::util::tokenize(option[1], stk::io::util::recognize(","));
00204             assert(tokens.size() == 3);
00205             offX = std::strtod(tokens[0].c_str(), NULL);
00206             offY = std::strtod(tokens[1].c_str(), NULL);
00207             offZ = std::strtod(tokens[2].c_str(), NULL);
00208           }
00209 
00210           else if (option[0] == "zdecomp") {
00211             // Option of the form  "zdecomp:1,1,2,2,1,2,...
00212             // Specifies the number of intervals in the z direction
00213             // for each processor.  The number of tokens must match
00214             // the number of processors.  Note that the new numZ will
00215             // be the sum of the intervals specified in this command.
00216             std::vector<std::string> tokens = stk::io::util::tokenize(option[1], stk::io::util::recognize(","));
00217             assert(tokens.size() == processorCount);
00218             std::vector<size_t> Zs;
00219             numZ = 0;
00220             for (size_t j = 0; j < processorCount; j++) {
00221               Zs.push_back(std::strtol(tokens[j].c_str(), NULL, 10));
00222               numZ += Zs[j];
00223             }
00224             myNumZ = Zs[myProcessor];
00225             myStartZ = 0;
00226             for (size_t j=0; j < myProcessor; j++) {
00227               myStartZ += Zs[j];
00228             }
00229           }
00230 
00231           else if (option[0] == "bbox") {
00232             // Bounding-Box Option of the form  "bbox:xmin,ymin,zmin,xmax,ymax,zmaxo
00233             std::vector<std::string> tokens = stk::io::util::tokenize(option[1], stk::io::util::recognize(","));
00234             assert(tokens.size() == 6);
00235             double xmin = std::strtod(tokens[0].c_str(), NULL);
00236             double ymin = std::strtod(tokens[1].c_str(), NULL);
00237             double zmin = std::strtod(tokens[2].c_str(), NULL);
00238             double xmax = std::strtod(tokens[3].c_str(), NULL);
00239             double ymax = std::strtod(tokens[4].c_str(), NULL);
00240             double zmax = std::strtod(tokens[5].c_str(), NULL);
00241 
00242             set_bbox(xmin, ymin, zmin,  xmax, ymax, zmax);
00243           }
00244 
00245           else if (option[0] == "rotate") {
00246             // Rotate Option of the form  "rotate:axis,angle,axis,angle,...
00247             std::vector<std::string> tokens = stk::io::util::tokenize(option[1], stk::io::util::recognize(","));
00248             assert(tokens.size() %2 == 0);
00249             for (size_t ir=0; ir < tokens.size();) {
00250               std::string axis = tokens[ir++];
00251               double angle_degree = std::strtod(tokens[ir++].c_str(), NULL);
00252               set_rotation(axis, angle_degree);
00253             }
00254           }
00255 
00256           else if (option[0] == "help") {
00257             std::cerr << "\nValid Options for GeneratedMesh parameter string:\n"
00258                       << "\tIxJxK -- specifies intervals; first option\n"
00259                       << "\toffset:xoff, yoff, zoff\n"
00260                       << "\tscale: xscl, yscl, zscl\n"
00261                       << "\tzdecomp:n1,n2,n3,...,n#proc\n"
00262                       << "\tbbox: xmin, ymin, zmin, xmax, ymax, zmax\n"
00263                       << "\trotate: axis,angle,axis,angle,...\n"
00264                       << "\tshell:xXyYzZ\n"
00265                       << "\tshow -- show mesh parameters\n"
00266                       << "\thelp -- show this list\n\n";
00267           }
00268 
00269           else if (option[0] == "show") {
00270       show_parameters();
00271           }
00272 
00273           else {
00274             std::cerr << "ERROR: Unrecognized option '" << option[0]
00275                       << "'.  It will be ignored.\n";
00276           }
00277         }
00278       }
00279 
00280       void GeneratedMesh::show_parameters() const
00281       {
00282   if (myProcessor == 0) {
00283     std::cerr << "\nMesh Parameters:\n"
00284         << "\tIntervals: " << numX << " by " << numY << " by " << numZ << "\n"
00285         << "\tX = " << sclX << " * (0.." << numX << ") + " << offX
00286         << "\tRange: " << offX << " <= X <= "  << offX + numX * sclX << "\n"
00287         << "\tY = " << sclY << " * (0.." << numY << ") + " << offY
00288         << "\tRange: " << offY << " <= Y <= "  << offY + numY * sclY << "\n"
00289         << "\tZ = " << sclZ << " * (0.." << numZ << ") + " << offZ
00290         << "\tRange: " << offZ << " <= Z <= "  << offZ + numZ * sclZ << "\n\n"
00291         << "\tNode Count (total)    = " << std::setw(9) << node_count() << "\n"
00292         << "\tElement Count (total) = " << std::setw(9) << element_count() << "\n"
00293         << "\tBlock Count           = " << std::setw(9) << block_count() << "\n\n";
00294     if (doRotation) {
00295       std::cerr << "\tRotation Matrix: \n\t" << std::scientific ;
00296       for (int ii=0; ii < 3; ii++) {
00297         for (int jj=0; jj < 3; jj++) {
00298     std::cerr << std::setw(14) << rotmat[ii][jj] << "\t";
00299         }
00300         std::cerr << "\n\t";
00301       }
00302       std::cerr << std::fixed << "\n";
00303     }
00304   }
00305       }
00306 
00307       size_t GeneratedMesh::node_count() const
00308       {
00309         return (numX+1) * (numY+1) * (numZ+1);
00310       }
00311 
00312       size_t GeneratedMesh::node_count_proc() const
00313       {
00314         return (numX+1) * (numY+1) * (myNumZ+1);
00315       }
00316 
00317       size_t GeneratedMesh::block_count() const
00318       {
00319         return shellBlocks.size() + 1;
00320       }
00321 
00322       size_t GeneratedMesh::element_count() const
00323       {
00324         size_t count = element_count(1);
00325         for (size_t i=0; i < shellBlocks.size(); i++) {
00326           count += element_count(i+2);
00327         }
00328         return count;
00329       }
00330 
00331       size_t GeneratedMesh::element_count_proc() const
00332       {
00333         size_t count = 0;
00334         for (size_t i=0; i < block_count(); i++) {
00335           count += element_count_proc(i+1);
00336         }
00337         return count;
00338       }
00339 
00340       size_t GeneratedMesh::element_count(size_t block_number) const
00341       {
00342         assert(block_number <= block_count());
00343 
00344         if (block_number == 1) {
00345           return numX * numY * numZ;
00346         } else {
00347           ShellLocation loc = shellBlocks[block_number-2];
00348           return shell_element_count(loc);
00349         }
00350       }
00351 
00352       size_t GeneratedMesh::shell_element_count(ShellLocation loc) const
00353       {
00354         switch (loc) {
00355           case MX:
00356           case PX:
00357             return numY * numZ;
00358           case MY:
00359           case PY:
00360             return numX * numZ;
00361           case MZ:
00362           case PZ:
00363             return numX * numY;
00364         }
00365         return 0;
00366       }
00367 
00368       size_t GeneratedMesh::element_count_proc(size_t block_number) const
00369       {
00370         assert(block_number <= block_count());
00371 
00372         if (block_number == 1) {
00373           return numX * numY * myNumZ;
00374         } else {
00375           ShellLocation loc = shellBlocks[block_number-2];
00376           return shell_element_count_proc(loc);
00377         }
00378       }
00379 
00380       size_t GeneratedMesh::shell_element_count_proc(ShellLocation loc) const
00381       {
00382         switch (loc) {
00383           case MX:
00384           case PX:
00385             return numY * myNumZ;
00386           case MY:
00387           case PY:
00388             return numX * myNumZ;
00389           case MZ:
00390             if (myProcessor == 0)
00391               return numX * numY;
00392             else
00393               return 0;
00394           case PZ:
00395             if (myProcessor == processorCount -1)
00396               return numX * numY;
00397             else
00398               return 0;
00399         }
00400         return 0;
00401       }
00402 
00403       std::pair<std::string,int> GeneratedMesh::topology_type(size_t block_number) const
00404       {
00405         assert(block_number <= block_count() && block_number > 0);
00406 
00407         if (block_number == 1) {
00408           return std::make_pair(std::string("hex8"), 8);;
00409         } else {
00410           return std::make_pair(std::string("shell4"), 4);
00411         }
00412       }
00413 
00414       void GeneratedMesh::node_map(std::vector<int> &map)
00415       {
00416         size_t count = node_count_proc();
00417         map.resize(count);
00418         size_t offset = myStartZ * (numX+1) * (numY+1);
00419         for (size_t i=0; i < count; i++) {
00420           map[i] = static_cast<int>(offset + i + 1);
00421         }
00422       }
00423 
00424       void GeneratedMesh::node_communication_map(std::vector<int> &map, std::vector<int> &proc)
00425       {
00426         size_t count = (numX+1) * (numY+1);
00427         size_t slab = count;
00428         if (myProcessor != 0 && myProcessor != processorCount-1)
00429           count *= 2;
00430 
00431         map.resize(count);
00432         proc.resize(count);
00433         size_t j = 0;
00434         if (myProcessor != 0) {
00435           size_t offset = myStartZ * (numX+1) * (numY+1);
00436           for (size_t i=0; i < slab; i++) {
00437             map[j] = static_cast<int>(offset + i + 1);
00438             proc[j++] = static_cast<int>(myProcessor-1);
00439           }
00440         }
00441         if (myProcessor != processorCount-1) {
00442           size_t offset = (myStartZ + myNumZ) * (numX+1) * (numY+1);
00443           for (size_t i=0; i < slab; i++) {
00444             map[j] = static_cast<int>(offset + i + 1);
00445             proc[j++] = static_cast<int>(myProcessor+1);
00446           }
00447         }
00448       }
00449 
00450       void GeneratedMesh::element_map(size_t block_number, std::vector<int> &map)
00451       {
00452         assert(block_number <= block_count() && block_number > 0);
00453 
00454         size_t count = element_count_proc(block_number);
00455         map.resize(count);
00456 
00457         if (block_number == 1) {
00458           // Hex block...
00459           count = element_count_proc(1);
00460           size_t offset = myStartZ * numX * numY;
00461           for (size_t i=0; i < count; i++) {
00462             map[i] = static_cast<int>(offset + i + 1);
00463           }
00464         } else {
00465           size_t start = element_count(1);
00466 
00467           // Shell blocks...
00468           for (size_t ib=0; ib < shellBlocks.size(); ib++) {
00469             count = element_count_proc(ib+2);
00470             if (static_cast<size_t>(block_number) == ib + 2) {
00471               size_t offset = 0;
00472               ShellLocation loc = shellBlocks[ib];
00473 
00474               switch (loc) {
00475               case MX:
00476               case PX:
00477                 offset = myStartZ * numY;
00478                 break;
00479 
00480               case MY:
00481               case PY:
00482                 offset = myStartZ * numX;
00483                 break;
00484 
00485               case MZ:
00486               case PZ:
00487                 offset = 0;
00488                 break;
00489               }
00490               for (size_t i=0; i < count; i++) {
00491                 map[i] = static_cast<int>(start + offset + i + 1);
00492               }
00493             } else {
00494               start += element_count(ib+2);
00495             }
00496           }
00497         }
00498       }
00499 
00500       void GeneratedMesh::element_map(std::vector<int> &map)
00501       {
00502         size_t count = element_count_proc();
00503         map.resize(count);
00504 
00505         size_t k = 0;
00506         // Hex block...
00507         count = element_count_proc(1);
00508         size_t offset = myStartZ * numX * numY;
00509         for (size_t i=0; i < count; i++) {
00510           map[k++] = static_cast<int>(offset + i + 1);
00511         }
00512 
00513         size_t start = element_count(1);
00514 
00515         // Shell blocks...
00516         for (size_t ib=0; ib < shellBlocks.size(); ib++) {
00517           count = element_count_proc(ib+2);
00518           offset = 0;
00519           ShellLocation loc = shellBlocks[ib];
00520 
00521           switch (loc) {
00522           case MX:
00523           case PX:
00524             offset = myStartZ * numY;
00525             break;
00526 
00527           case MY:
00528           case PY:
00529             offset = myStartZ * numX;
00530             break;
00531 
00532           case MZ:
00533           case PZ:
00534             offset = 0;
00535             break;
00536           }
00537           for (size_t i=0; i < count; i++) {
00538             map[k++] = static_cast<int>(start + offset + i + 1);
00539           }
00540           start += element_count(ib+2);
00541         }
00542       }
00543 
00544       void GeneratedMesh::element_surface_map(ShellLocation loc, std::vector<int> &map)
00545       {
00546         size_t count = shell_element_count_proc(loc);
00547         map.resize(2*count);
00548 
00549         size_t offset = 0;
00550 
00551         switch (loc) {
00552         case MX:
00553           for( size_t k = 0, index = 0; k < numZ; ++k )
00554             for( size_t j = 0; j < numY; ++j, ++index )
00555             {
00556               size_t ielem = k*numX*numY+j*numX;
00557               map[2*index] = ielem + 1; // 1-based elem id
00558               map[2*index+1] = 3; // 0-based local face id
00559             }
00560           break;
00561 
00562         case PX:
00563           offset = numX-1;
00564           for( size_t k = 0, index = 0; k < numZ; ++k )
00565             for( size_t j = 0; j < numY; ++j, ++index )
00566             {
00567               size_t ielem = offset + k*numX*numY+j*numX + 1;
00568               map[2*index] = ielem;
00569               map[2*index+1] = 1; // 0-based local face id
00570             }
00571           break;
00572 
00573         case MY:
00574           for( size_t k = 0, index = 0; k < numZ; ++k )
00575             for( size_t i = 0; i < numX; ++i, ++index )
00576             {
00577               size_t ielem = k*numX*numY+i + 1;
00578               map[2*index] = ielem;
00579               map[2*index+1] = 0; // 0-based local face id
00580             }
00581           break;
00582         case PY:
00583           offset = (numY-1)*numX;
00584           for( size_t k = 0, index = 0; k < numZ; ++k )
00585             for( size_t i = 0; i < numX; ++i, ++index )
00586             {
00587               size_t ielem = offset + k*numX*numY+i;
00588               map[2*index] = ielem + 1;
00589               map[2*index+1] = 2; // 0-based local face id
00590             }
00591           break;
00592 
00593         case MZ:
00594           for( size_t j = 0; j < numY; ++j )
00595             for( size_t i = 0; i < numX; ++i )
00596             {
00597               size_t ielem = j*numX+i;
00598               size_t index = 2*ielem;
00599               map[index] = ielem + 1;
00600               map[index+1] = 4; // 0-based local face id
00601             }
00602           break;
00603 
00604         case PZ:
00605           offset = (numZ-1)*numX*numY;
00606           for( size_t j = 0; j < numY; ++j )
00607             for( size_t i = 0; i < numX; ++i )
00608             {
00609               size_t ielem = offset + j*numX+i;
00610               size_t index = 2*(ielem-offset);
00611               map[index] = ielem + 1;
00612               map[index+1] = 5; // 0-based local face id
00613             }
00614           break;
00615         }
00616       }
00617 
00618       void GeneratedMesh::coordinates(std::vector<double> &coord) const
00619       {
00620         /* create global coordinates */
00621         size_t count = node_count_proc();
00622         coord.resize(count * 3);
00623 
00624         size_t k = 0;
00625         for (size_t m=myStartZ; m < myStartZ+myNumZ+1; m++) {
00626           for (size_t i=0; i < numY+1; i++) {
00627             for (size_t j=0; j < numX+1; j++) {
00628               coord[k++] = sclX * static_cast<double>(j) + offX;
00629               coord[k++] = sclY * static_cast<double>(i) + offY;
00630               coord[k++] = sclZ * static_cast<double>(m) + offZ;
00631             }
00632           }
00633         }
00634 
00635         if (doRotation) {
00636           for (size_t i=0; i < count*3; i+=3) {
00637             double xn = coord[i+0];
00638             double yn = coord[i+1];
00639             double zn = coord[i+2];
00640             coord[i+0] = xn * rotmat[0][0] + yn * rotmat[1][0] + zn * rotmat[2][0];
00641             coord[i+1] = xn * rotmat[0][1] + yn * rotmat[1][1] + zn * rotmat[2][1];
00642             coord[i+2] = xn * rotmat[0][2] + yn * rotmat[1][2] + zn * rotmat[2][2];
00643           }
00644         }
00645       }
00646 
00647       void GeneratedMesh::coordinates(std::vector<double> &x,
00648                                       std::vector<double> &y,
00649                                       std::vector<double> &z) const
00650       {
00651         /* create global coordinates */
00652         size_t count = node_count_proc();
00653         x.resize(count);
00654         y.resize(count);
00655         z.resize(count);
00656 
00657         size_t k = 0;
00658         for (size_t m=myStartZ; m < myStartZ+myNumZ+1; m++) {
00659           for (size_t i=0; i < numY+1; i++) {
00660             for (size_t j=0; j < numX+1; j++) {
00661               x[k] = sclX * static_cast<double>(j) + offX;
00662               y[k] = sclY * static_cast<double>(i) + offY;
00663               z[k] = sclZ * static_cast<double>(m) + offZ;
00664               ++k;
00665             }
00666           }
00667         }
00668         if (doRotation) {
00669           for (size_t i=0; i < count; i++) {
00670             double xn = x[i];
00671             double yn = y[i];
00672             double zn = z[i];
00673             x[i] = xn * rotmat[0][0] + yn * rotmat[1][0] + zn * rotmat[2][0];
00674             y[i] = xn * rotmat[0][1] + yn * rotmat[1][1] + zn * rotmat[2][1];
00675             z[i] = xn * rotmat[0][2] + yn * rotmat[1][2] + zn * rotmat[2][2];
00676           }
00677         }
00678       }
00679 
00680       void GeneratedMesh::connectivity(size_t block_number, std::vector<int> &connect) const
00681       {
00682         assert(block_number <= block_count());
00683 
00684         size_t xp1yp1 = (numX+1) * (numY+1);
00685 
00686         /* build connectivity array (node list) for mesh */
00687         if (block_number == 1) {  // HEX Element Block
00688           connect.resize(element_count_proc(block_number)*8);
00689 
00690           size_t cnt = 0;
00691           for (size_t m=myStartZ; m < myNumZ+myStartZ; m++) {
00692             for (size_t i=0, k=0; i < numY; i++) {
00693               for (size_t j=0; j < numX; j++, k++) {
00694                 size_t base = (m*xp1yp1) + k + i + 1 - myStartZ * xp1yp1;
00695                 ;
00696                 connect[cnt++] = static_cast<int>(base);
00697                 connect[cnt++] = static_cast<int>(base+1);
00698                 connect[cnt++] = static_cast<int>(base+numX+2);
00699                 connect[cnt++] = static_cast<int>(base+numX+1);
00700 
00701                 connect[cnt++] = static_cast<int>(xp1yp1 + base);
00702                 connect[cnt++] = static_cast<int>(xp1yp1 + base+1);
00703                 connect[cnt++] = static_cast<int>(xp1yp1 + base+numX+2);
00704                 connect[cnt++] = static_cast<int>(xp1yp1 + base+numX+1);
00705               }
00706             }
00707           }
00708         } else { // Shell blocks....
00709           ShellLocation loc = shellBlocks[block_number-2];
00710           connect.resize(element_count(block_number)*4);
00711 
00712           size_t cnt = 0;
00713           switch (loc) {
00714           case MX:  // Minumum X Face
00715             for (size_t i=0; i < myNumZ; i++) {
00716               size_t layer_off = i * xp1yp1;
00717               for (size_t j=0; j < numY; j++) {
00718                 size_t base = layer_off + j * (numX+1) + 1;
00719                 connect[cnt++] = static_cast<int>(base);
00720                 connect[cnt++] = static_cast<int>(base + xp1yp1);
00721                 connect[cnt++] = static_cast<int>(base + xp1yp1 + (numX+1));
00722                 connect[cnt++] = static_cast<int>(base + (numX+1));
00723               }
00724             }
00725             break;
00726           case PX: // Maximum X Face
00727             for (size_t i=0; i < myNumZ; i++) {
00728               size_t layer_off = i * xp1yp1;
00729               for (size_t j=0; j < numY; j++) {
00730                 size_t base = layer_off + j * (numX+1) + numX + 1;
00731                 connect[cnt++] = static_cast<int>(base);
00732                 connect[cnt++] = static_cast<int>(base + (numX+1));
00733                 connect[cnt++] = static_cast<int>(base + xp1yp1 + (numX+1));
00734                 connect[cnt++] = static_cast<int>(base + xp1yp1);
00735               }
00736             }
00737             break;
00738           case MY: // Minumum Y Face
00739             for (size_t i=0; i < myNumZ; i++) {
00740               size_t layer_off = i * xp1yp1;
00741               for (size_t j=0; j < numX; j++) {
00742                 size_t base = layer_off + j + 1;
00743                 connect[cnt++] = static_cast<int>(base);
00744                 connect[cnt++] = static_cast<int>(base + 1);
00745                 connect[cnt++] = static_cast<int>(base + xp1yp1 + 1);
00746                 connect[cnt++] = static_cast<int>(base + xp1yp1);
00747               }
00748             }
00749             break;
00750           case PY: // Maximum Y Face
00751             for (size_t i=0; i < myNumZ; i++) {
00752               size_t layer_off = i * xp1yp1;
00753               for (size_t j=0; j < numX; j++) {
00754                 size_t base = layer_off + (numX+1)*(numY) + j + 1;
00755                 connect[cnt++] = static_cast<int>(base);
00756                 connect[cnt++] = static_cast<int>(base + xp1yp1);
00757                 connect[cnt++] = static_cast<int>(base + xp1yp1 + 1);
00758                 connect[cnt++] = static_cast<int>(base + 1);
00759               }
00760             }
00761             break;
00762           case MZ: // Minumum Z Face
00763             if (myProcessor == 0) {
00764               for (size_t i=0, k=0; i < numY; i++) {
00765                 for (size_t j=0; j < numX; j++, k++) {
00766                   size_t base = i + k + 1;
00767                   connect[cnt++] = static_cast<int>(base);
00768                   connect[cnt++] = static_cast<int>(base+numX+1);
00769                   connect[cnt++] = static_cast<int>(base+numX+2);
00770                   connect[cnt++] = static_cast<int>(base+1);
00771                 }
00772               }
00773             }
00774             break;
00775           case PZ: // Maximum Z Face
00776             if (myProcessor == processorCount-1) {
00777               for (size_t i=0, k=0; i < numY; i++) {
00778                 for (size_t j=0; j < numX; j++, k++) {
00779                   size_t base = xp1yp1 * (numZ - myStartZ) + k + i + 1;
00780                   connect[cnt++] = static_cast<int>(base);
00781                   connect[cnt++] = static_cast<int>(base+1);
00782                   connect[cnt++] = static_cast<int>(base+numX+2);
00783                   connect[cnt++] = static_cast<int>(base+numX+1);
00784                 }
00785               }
00786             }
00787             break;
00788           }
00789         }
00790         return;
00791       }
00792 
00793       void GeneratedMesh::set_rotation(const std::string &axis, double angle_degrees)
00794       {
00795         // PI / 180. Used in converting angle in degrees to radians
00796         static double degang = std::atan2(0.0, -1.0) / 180.0;
00797 
00798         doRotation = true;
00799 
00800         int n1 = -1;
00801   int n2 = -1;
00802   int n3 = -1;
00803 
00804         if (axis == "x" || axis == "X") {
00805           n1 = 1; n2 = 2; n3 = 0;
00806         } else if (axis == "y" || axis == "Y") {
00807           n1 = 2; n2 = 0; n3 = 1;
00808         } else if (axis == "z" || axis == "Z") {
00809           n1 = 0; n2 = 1; n3 = 2;
00810         } else {
00811           std::cerr << "\nInvalid axis specification '" << axis << "'. Valid options are 'x', 'y', or 'z'\n";
00812           return;
00813         }
00814 
00815         double ang = angle_degrees * degang; // Convert angle in degrees to radians
00816         double cosang = std::cos(ang);
00817         double sinang = std::sin(ang);
00818 
00819   assert(n1 >= 0 && n2 >= 0 && n3 >= 0);
00820         double by[3][3];
00821         by[n1][n1] =  cosang;
00822         by[n2][n1] = -sinang;
00823         by[n1][n3] = 0.0;
00824         by[n1][n2] =  sinang;
00825         by[n2][n2] =  cosang;
00826         by[n2][n3] = 0.0;
00827         by[n3][n1] = 0.0;
00828         by[n3][n2] = 0.0;
00829         by[n3][n3] = 1.0;
00830 
00831         double res[3][3];
00832         for (size_t i=0; i < 3; i++) {
00833           res[i][0] = rotmat[i][0]*by[0][0] + rotmat[i][1]*by[1][0] + rotmat[i][2]*by[2][0];
00834           res[i][1] = rotmat[i][0]*by[0][1] + rotmat[i][1]*by[1][1] + rotmat[i][2]*by[2][1];
00835           res[i][2] = rotmat[i][0]*by[0][2] + rotmat[i][1]*by[1][2] + rotmat[i][2]*by[2][2];
00836         }
00837 
00838 #if 1
00839         std::memcpy(rotmat, res, 9*sizeof(double));
00840 #else
00841         for (int i=0; i < 3; i++) {
00842           for (int j=0; j < 3; j++) {
00843             rotmat[i][j] = res[i][j];
00844           }
00845         }
00846 #endif
00847       }
00848     }
00849   }
00850 }
00851 
00852 #if defined(DEBUG)
00853 #include <sstream>
00854   std::string ToString(size_t t) {
00855     std::ostringstream os;
00856     os << t;
00857     return os.str();
00858   }
00859 
00860 #include <exodusII.h>
00861 #include <ne_nemesisI.h>
00862 
00863   int main() {
00864     int num_processors = 8;
00865     for (int proc = 0; proc < num_processors; proc++) {
00866 
00867       stk::GeneratedMesh mesh(100, 125, 10*num_processors, num_processors, proc);
00868 
00869       std::cerr << "Node Count (total)    = " << mesh.node_count() << "\n";
00870       std::cerr << "Node Count (proc)     = " << mesh.node_count_proc() << "\n";
00871       std::cerr << "Element Count (total) = " << mesh.element_count() << "\n";
00872       std::cerr << "Element Count (proc)  = " << mesh.element_count_proc() << "\n";
00873       std::cerr << "Block Count           = " << mesh.block_count() << "\n";
00874 
00875       int CPU_word_size = 8;                   /* sizeof(float) */
00876       int IO_word_size = 8;                    /* (4 bytes) */
00877       std::string name = "test-scale.e";
00878       if (num_processors > 1) {
00879         name += "." + ToString(num_processors) + "." + ToString(proc);
00880       }
00881       int exoid = ex_create (name.c_str(), EX_CLOBBER, &CPU_word_size, &IO_word_size);
00882 
00883       int num_nodes = mesh.node_count_proc();
00884       int num_elems = mesh.element_count_proc();
00885       int num_elem_blk = mesh.block_count();
00886       int error = ex_put_init (exoid, "title", 3, num_nodes, num_elems,
00887                                num_elem_blk, 0, 0);
00888 
00889       if (num_processors > 1) {
00890         std::vector<int> nodes;
00891         std::vector<int> procs;
00892         mesh.node_communication_map(nodes, procs);
00893 
00894         int node_map_ids[1] = {1};
00895         int node_map_node_cnts[1] = {procs.size()};
00896         ne_put_init_info(exoid, num_processors, 1, "p");
00897         ne_put_loadbal_param(exoid, 0, 0, 0, 0, 0, 1, 0, proc);
00898         ne_put_cmap_params(exoid, node_map_ids, node_map_node_cnts, 0, 0, proc);
00899         ne_put_node_cmap(exoid, 1, &nodes[0], &procs[0], proc);
00900       }
00901 
00902       for (int i=1; i < mesh.block_count(); i++) {
00903         std::cerr << "Block " << i+1 << " has " << mesh.element_count_proc(i+1) << " "
00904                   << mesh.topology_type(i+1).first <<" elements\n";
00905 
00906         std::string btype = mesh.topology_type(i+1).first;
00907         error = ex_put_elem_block(exoid, i+1, btype.c_str(),
00908                                   mesh.element_count_proc(i+1),
00909                                   mesh.topology_type(i+1).second, 0);
00910       }
00911       {
00912         std::cerr << "Block " << 1 << " has " << mesh.element_count_proc(1) << " "
00913                   << mesh.topology_type(1).first <<" elements\n";
00914 
00915         std::string btype = mesh.topology_type(1).first;
00916         error = ex_put_elem_block(exoid, 1, btype.c_str(),
00917                                   mesh.element_count_proc(1),
00918                                   mesh.topology_type(1).second, 0);
00919       }
00920 
00921       if (num_processors > 1) {
00922         std::vector<int> map;
00923         mesh.node_map(map);
00924         ex_put_id_map(exoid, EX_NODE_MAP, &map[0]);
00925 
00926         mesh.element_map(map);
00927         ex_put_id_map(exoid, EX_ELEM_MAP, &map[0]);
00928       }
00929 
00930       std::cerr << "Outputting connectivity...\n";
00931       for (int i=1; i < mesh.block_count(); i++) {
00932         if (mesh.element_count_proc(i+1) > 0) {
00933           std::vector<int> connectivity;
00934           mesh.connectivity(i+1, connectivity);
00935           ex_put_elem_conn(exoid, i+1, &connectivity[0]);
00936         }
00937       }
00938       {
00939         std::vector<int> connectivity;
00940         mesh.connectivity(1, connectivity);
00941         ex_put_elem_conn(exoid, 1, &connectivity[0]);
00942       }
00943 
00944       {
00945         std::vector<double> x;
00946         std::vector<double> y;
00947         std::vector<double> z;
00948         mesh.coordinates(x, y, z);
00949         error = ex_put_coord (exoid, &x[0], &y[0], &z[0]);
00950       }
00951 
00952       ex_close(exoid);
00953     }
00954   }
00955 #endif
Generated on Wed Apr 13 10:05:48 2011 for Sierra Toolkit by  doxygen 1.6.3