MeshKit  1.0
AssyGen.cpp
Go to the documentation of this file.
00001 #include "meshkit/AssyGen.hpp"
00002 
00003 namespace MeshKit
00004 {
00005   // static registration of this  mesh scheme
00006   moab::EntityType AssyGen_tps[] = { moab::MBVERTEX,
00007                                      moab::MBEDGE,
00008                                      moab::MBTRI,
00009                                      moab::MBHEX,
00010                                      moab::MBMAXTYPE};
00011   const moab::EntityType* AssyGen::output_types()
00012   { return AssyGen_tps; }
00013 
00014   AssyGen::AssyGen( MKCore *mk, const MEntVector &me_vec)
00015     : MeshScheme( mk, me_vec),
00016       igeomImpl(mk->igeom_instance())
00017   {
00018     m_nPlanar = 0; //default is 3D
00019     m_nLineNumber = 0;
00020     root_set= NULL;
00021     szComment = "!";
00022     MAXCHARS = 300;
00023     pi = M_PI;
00024     m_dRadialSize = -1.0;
00025     m_dAxialSize = -1.0;
00026     m_dTetMeshSize = -1.0;
00027     m_nDimensions = 0;
00028     m_nMaterialSetId = 1;
00029     m_nNeumannSetId = 1;
00030     m_szEngine = "acis";
00031     m_szMeshType = "hex";
00032     m_nDuct = 0;
00033     m_nDuctNum = 0;
00034     m_nJouFlag = 0;
00035     m_szSideset = "yes";
00036     m_dMergeTol = 1e-4;
00037   }
00038 
00039   AssyGen::~AssyGen()
00040   {}
00041 
00042   bool AssyGen::add_modelent(ModelEnt *model_ent)
00043   {
00044     return MeshOp::add_modelent(model_ent);
00045   }
00046 
00047   void AssyGen::setup_this()
00048   {
00049 
00050     // start the timer 
00051     CClock Timer;
00052     clock_t sTime = clock();
00053     std::string szDateTime;
00054     Timer.GetDateTime (szDateTime);
00055     std::cout << "\nStarting out at : " << szDateTime << "\n";
00056   
00057     //count pin cylinders and cell material, needed for setting array size before actual read
00058     ReadInputPhase1 ();
00059 
00060     // read the problem size and create pincell
00061     ReadAndCreate ();
00062 
00063     // create the .jou file
00064     CreateCubitJournal();
00065 
00066     // get the current date and time
00067     Timer.GetDateTime (szDateTime);
00068     std::cout << "Ending at : " << szDateTime;
00069  
00070     // compute the elapsed time
00071     std::cout << "Elapsed wall clock time: " << Timer.DiffTime ()
00072               << " seconds or " << (Timer.DiffTime ())/60.0 << " mins\n";
00073 
00074     std::cout << "Total CPU time used: " << (double) (clock() - sTime)/CLOCKS_PER_SEC \
00075               << " seconds" << std::endl; 
00076   }
00077 
00078   void AssyGen::execute_this()
00079   {
00080   }
00081 
00082   void AssyGen::PrepareIO (int argc, char *argv[], std::string  TestDir)
00083   // ---------------------------------------------------------------------------
00084   // Function: Obtains file names and opens input/output files
00085   // Input:    command line arguments
00086   // Output:   none
00087   // ---------------------------------------------------------------------------
00088   {
00089     std::cout << '\n';
00090     std::cout << "\t\t---------------------------------------------------------" << '\n';
00091     std::cout << "\t\tProgram to Generate Nuclear Reactor Assembly Geometries      " << '\n';
00092     std::cout << "\t\t\t\tArgonne National Laboratory" << '\n';
00093     std::cout << "\t\t\t\t        2010-2014         " << '\n';
00094     std::cout << "\t\t---------------------------------------------------------" << '\n';
00095     std::cout << "\nsee README file for using the program and details on various cards.\n"<< std::endl;
00096 
00097     // set and open input output files
00098     bool bDone = false;
00099     do{
00100       if (2 == argc) {
00101         m_szFile = argv[1];
00102         m_szInFile=m_szFile+".inp";
00103         m_szJouFile = m_szFile+".jou";
00104         m_szSchFile = m_szFile+".template.jou";
00105       }
00106       else if (3 == argc) {
00107         int i=1;// will loop through arguments, and process them
00108         for (i=1; i<argc-1 ; i++) {
00109           if (argv[i][0]=='-') {
00110             switch (argv[i][1])
00111               {
00112               case 'j':
00113                 {
00114                   m_nJouFlag = 1;
00115                   std::cout << "Creating journal file only.\n Geometry file must exist in the same directory." << std::endl;
00116                   m_szFile = argv[2];
00117                   m_szInFile=m_szFile+".inp";
00118                   m_szJouFile = m_szFile+".jou";
00119                   m_szSchFile = m_szFile+".template.jou";
00120                   break;
00121                 }
00122               case 'h':
00123                 {
00124                   std::cout << "\nInstruction on writing assygen input file can also be found at: " << std::endl;
00125                   std::cout << "        https://trac.mcs.anl.gov/projects/fathom/browser/MeshKit/trunk/rgg/README" << std::endl;
00126                   std::cout << "Usage: assygen [-j -h] <input file name without extension>"<< std::endl;
00127                   std::cout << "        -j create journal file only" << std::endl;
00128                   std::cout << "        -h print help" << std::endl;
00129 
00130                   exit(0);
00131                   break;
00132                 }
00133               }
00134           }
00135         }
00136       }
00137       else if (1 == argc){
00138         std::cout << "\nInstruction on writing assygen input file can also be found at: " << std::endl;
00139         std::cout << "        https://trac.mcs.anl.gov/projects/fathom/browser/MeshKit/trunk/rgg/README" << std::endl;
00140         std::cout << "Usage: assygen [-t -j -h] <input file name without extension>"<< std::endl;
00141         std::cout << "        -t print timing and memory usage info in each step" << std::endl;
00142         std::cout << "        -j create journal file only" << std::endl;
00143         std::cout << "        -h print help" << std::endl;
00144         std::cout << "\nRunning default case:\n" << std::endl;
00145 
00146         m_szInFile =  TestDir + "/" + DEFAULT_TEST_FILE;
00147         m_szGeomFile = (char *)TEST_FILE_NAME;
00148         m_szJouFile = (char *)TEST_FILE_NAME;
00149         m_szFile =  (char *)TEST_FILE_NAME;
00150         m_szInFile+=".inp";
00151         m_szJouFile+=".jou";
00152         m_szSchFile = m_szFile+".template.jou";
00153         std::cout <<"  No file specified.  Defaulting to: " << m_szInFile
00154                   << "  " << m_szJouFile << std::endl;
00155       }
00156       // open the file
00157       m_FileInput.open (m_szInFile.c_str(), std::ios::in);
00158       if (!m_FileInput){
00159         std::cout << "Unable to open file" << std::endl;
00160         std::cout << "Usage: assygen <input filename WITHOUT EXTENSION>"<< std::endl;
00161         m_FileInput.clear ();
00162         exit(1);
00163       }
00164       else
00165         bDone = true; // file opened successfully
00166     } while (!bDone);
00167     std::cout << "\nEntered input file name: " <<  m_szInFile <<std::endl;
00168 
00169     // open the file
00170     do{
00171       m_FileOutput.open (m_szJouFile.c_str(), std::ios::out);
00172       if (!m_FileOutput){
00173         std::cout << "Unable to open o/p file" << std::endl;
00174         m_FileOutput.clear ();
00175         exit(1);
00176       }
00177       else
00178         bDone = true; // file opened successfully
00179     } while (!bDone);
00180 
00181     // open the template journal file for writing
00182     do{
00183       m_SchemesFile.open (m_szSchFile.c_str(), std::ios::out);
00184       if (!m_SchemesFile){
00185         std::cout << "Unable to open o/p file" << std::endl;
00186         m_SchemesFile.clear ();
00187         exit(1);
00188       }
00189       else
00190         bDone = true; // file opened successfully
00191     } while (!bDone);
00192 
00193     std::cout<<"\no/p Cubit journal file name: "<< m_szJouFile
00194              << std::endl;
00195 
00196   }
00197 
00198   void AssyGen::ReadInputPhase1 ()
00199   // -------------------------------------------------------------------------------------------
00200   // Function: reads the input file to count the no. of cyl in a pincell, before the actual read
00201   // Input:    none
00202   // Output:   none
00203   // -------------------------------------------------------------------------------------------
00204   {
00205     CParser Parse;
00206     int nCyl =0, nCellMat=0, nInputLines=0;
00207     std::string card, szVolId, szVolAlias;
00208     // count the total number of cylinder commands in each pincell
00209     for(;;){
00210       if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
00211                                MAXCHARS, szComment))
00212         IOErrorHandler (INVALIDINPUT);
00213 
00214       if (szInputString.substr(0,10) == "geomengine"){
00215         std::istringstream szFormatString (szInputString);
00216         szFormatString >> card >> m_szEngine;
00217         if( ((strcmp (m_szEngine.c_str(), "acis") != 0) &&
00218              (strcmp (m_szEngine.c_str(), "occ") != 0)) || szFormatString.fail())
00219           IOErrorHandler(EGEOMENGINE);
00220       }
00221       if (szInputString.substr(0,8) == "meshtype"){
00222         std::istringstream szFormatString (szInputString);
00223         szFormatString >> card >> m_szMeshType;
00224         if( ((strcmp (m_szMeshType.c_str(), "hex") != 0) &&
00225              (strcmp (m_szMeshType.c_str(), "tet") != 0)) || szFormatString.fail())
00226           IOErrorHandler(INVALIDINPUT);
00227       }
00228       if (szInputString.substr(0,4) == "duct" || szInputString.substr(0,10) == "dimensions"){
00229         ++m_nDuct;
00230       }
00231       if (szInputString.substr(0,8) == "pincells"){
00232         std::istringstream szFormatString (szInputString);
00233         szFormatString >> card >> m_nPincells;
00234         if(m_nPincells>0)
00235           m_Pincell.SetSize(m_nPincells);
00236         else if(m_nPincells ==0)
00237           m_Pincell.SetSize(1); // assume for using dummy pincell
00238 
00239         // count the number of cylinder lines for each pincell
00240         for (int i=1; i<=m_nPincells; i++){
00241           // read the no. of input lines first pincell
00242           if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
00243                                    MAXCHARS, szComment))
00244             IOErrorHandler (INVALIDINPUT);
00245           std::istringstream szFormatString1 (szInputString);
00246           szFormatString1 >> szVolId >> szVolAlias >> nInputLines;
00247           if(szFormatString1.fail())
00248             IOErrorHandler(INVALIDINPUT);
00249           // loop thru the input lines of each pincell
00250           for(int l=1; l<=nInputLines; l++){
00251             if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
00252                                      MAXCHARS, szComment))
00253               IOErrorHandler (INVALIDINPUT);
00254             if (szInputString.substr(0,8) == "cylinder"){
00255               ++nCyl;
00256             }
00257             if (szInputString.substr(0,12) == "cellmaterial"){
00258               ++nCellMat;
00259             }
00260           }
00261 
00262           // set the sizes
00263           if(nCyl>0){
00264             if  (nCellMat!=0){
00265               m_Pincell(i).SetCellMatSize(nCyl);
00266             }
00267             m_Pincell(i).SetNumCyl(nCyl);
00268           }
00269           else if(nCyl ==0){
00270             if(nInputLines >0)
00271               m_Pincell(i).SetCellMatSize(nCellMat);
00272           }
00273           nCyl = 0;
00274           nCellMat = 0;
00275         }
00276       }
00277       // breaking condition
00278       if(szInputString.substr(0,3) == "end"){
00279         std::istringstream szFormatString (szInputString);
00280         break;
00281       }
00282     }
00283 
00284     //ACIS ENGINE
00285 #ifdef HAVE_ACIS
00286     //  if(m_szEngine == "acis"){
00287     m_szGeomFile = m_szFile+".sat";
00288     //  }
00289 #elif defined(HAVE_OCC)
00290     //  OCC ENGINE
00291     //  if (m_szEngine == "occ"){
00292     m_szGeomFile = m_szFile+".stp";
00293     //  }
00294 #endif
00295     std::cout << "\no/p geometry file name: " <<  m_szGeomFile <<std::endl;
00296 
00297 
00298   }
00299 
00300   void AssyGen::ReadPinCellData ( int i)
00301   //---------------------------------------------------------------------------
00302   //Function: reading pincell i from file and storing the data
00303   //Input:    none
00304   //Output:   none
00305   //---------------------------------------------------------------------------
00306   {
00307     CParser Parse;
00308     std::string card, szVolId, szVolAlias, szIFlag;
00309     int nInputLines, nMaterials, nCyl = 0, nRadii=0, nCellMat=0;
00310     double dLZ=0.0, dFlatF=0.0, dPX=0.0, dPY=0.0, dPZ=0.0;
00311     CVector <std::string> szVMatName, szVMatAlias, szVCylMat, szVCellMat;
00312     CVector<double> dVCoor(2), dVCylRadii, dVCylZPos, dZVStart, dZVEnd;
00313 
00314     //loop over input lines
00315     if (m_szGeomType == "rectangular"){
00316 
00317       std::cout << "\ngetting volume id";
00318       if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
00319                                MAXCHARS, szComment))
00320         IOErrorHandler (INVALIDINPUT);
00321       std::istringstream szFormatString (szInputString);
00322       szFormatString >> szVolId >> szVolAlias >> nInputLines >> szIFlag;
00323 
00324       // error checking
00325       if( (strcmp (szVolAlias.c_str(), "") == 0) ||
00326           (strcmp (szVolId.c_str(), "") == 0))
00327         IOErrorHandler(EPIN);
00328       if( nInputLines < 0 )
00329         IOErrorHandler(ENEGATIVE);
00330 
00331       m_Pincell(i).SetLineOne (szVolId, szVolAlias, nInputLines);
00332       if(szIFlag == "intersect"){
00333         m_Pincell(i).SetIntersectFlag(1);
00334       }
00335       else{
00336         m_Pincell(i).SetIntersectFlag(0);
00337       }
00338       for(int l=1; l<=nInputLines; l++){
00339         if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
00340                                  MAXCHARS, szComment))
00341           IOErrorHandler (INVALIDINPUT);
00342         if (szInputString.substr(0,5) == "pitch"){
00343 
00344           std::istringstream szFormatString (szInputString);
00345           std::cout << "\ngetting pitch data";
00346           szFormatString >> card >> dPX >> dPY >> dPZ;
00347 
00348           if( dPX < 0 || dPY < 0 || dPZ < 0 || szFormatString.fail())
00349             IOErrorHandler(ENEGATIVE);
00350           m_Pincell(i).SetPitch (dPX, dPY, dPZ);
00351         }
00352         if (szInputString.substr(0,9) == "materials"){
00353 
00354           std::istringstream szFormatString (szInputString);
00355           szFormatString >> card >> nMaterials;
00356           if(szFormatString.fail())
00357             IOErrorHandler(INVALIDINPUT);
00358 
00359           //setting local arrays
00360           szVMatName.SetSize(nMaterials);
00361           szVMatAlias.SetSize(nMaterials);
00362 
00363           //set class variable sizes
00364           m_Pincell(i).SetMatArray(nMaterials);
00365           std::cout << "\ngetting material data";
00366           for(int j=1; j<= nMaterials; j++){
00367             szFormatString >> szVMatName(j) >> szVMatAlias(j);
00368             if(szFormatString.fail())
00369               IOErrorHandler(INVALIDINPUT);
00370           }
00371           m_Pincell(i).SetMat(szVMatName, szVMatAlias);
00372         }
00373         if (szInputString.substr(0,8) == "cylinder"){
00374 
00375           ++nCyl;
00376           std::cout << "\ngetting cylinder data";
00377           std::istringstream szFormatString (szInputString);
00378           szFormatString >> card >> nRadii >> dVCoor(1) >> dVCoor(2);
00379           if(szFormatString.fail())
00380             IOErrorHandler(INVALIDINPUT);
00381           m_Pincell(i).SetCylSizes(nCyl, nRadii);
00382           m_Pincell(i).SetCylPos(nCyl, dVCoor);
00383 
00384           //set local array
00385           dVCylRadii.SetSize(nRadii);
00386           szVCylMat.SetSize(nRadii);
00387           dVCylZPos.SetSize(2);
00388           m_Pincell(i).SetCylSizes(nCyl, nRadii);
00389 
00390           // reading ZCoords
00391           for(int k=1; k<=2; k++){
00392             szFormatString >> dVCylZPos(k);
00393             if(szFormatString.fail())
00394               IOErrorHandler(INVALIDINPUT);
00395           }
00396           m_Pincell(i).SetCylZPos(nCyl, dVCylZPos);
00397 
00398           // reading Radii
00399           for(int l=1; l<= nRadii; l++){
00400             szFormatString >> dVCylRadii(l);
00401             if( dVCylRadii(l) < 0  || szFormatString.fail())
00402               IOErrorHandler(ENEGATIVE);
00403           }
00404           m_Pincell(i).SetCylRadii(nCyl, dVCylRadii);
00405 
00406           // reading Material alias
00407           for(int m=1; m<= nRadii; m++){
00408             szFormatString >> szVCylMat(m);
00409             if(strcmp (szVCylMat(m).c_str(), "") == 0 || szFormatString.fail())
00410               IOErrorHandler(EALIAS);
00411           }
00412           m_Pincell(i).SetCylMat(nCyl, szVCylMat);
00413         }
00414         if (szInputString.substr(0,12) == "cellmaterial"){
00415 
00416           std::cout << "\ngetting cell material data\n";
00417           std::istringstream szFormatString (szInputString);
00418           szFormatString >> card;
00419 
00420           //set local arrays
00421           m_Pincell(i).GetCellMatSize(nCellMat); // since size of cell material already set equal to number of cylinders
00422           dZVStart.SetSize(nCellMat);
00423           dZVEnd.SetSize(nCellMat);
00424           szVCellMat.SetSize(nCellMat);
00425 
00426           for(int k=1; k<=nCellMat; k++){
00427             szFormatString >> dZVStart(k)>> dZVEnd(k) >> szVCellMat(k);
00428             if(strcmp (szVCellMat(k).c_str(), "") == 0 || szFormatString.fail())
00429               IOErrorHandler(EALIAS);
00430           }
00431           m_Pincell(i).SetCellMat(dZVStart, dZVEnd, szVCellMat);
00432         }
00433       }
00434     }//if rectangular ends
00435 
00436     if (m_szGeomType == "hexagonal"){
00437 
00438       std::cout << "\ngetting volume id";
00439       if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
00440                                MAXCHARS, szComment))
00441         IOErrorHandler (INVALIDINPUT);
00442       std::istringstream szFormatString (szInputString);
00443       szFormatString >> szVolId >> szVolAlias >> nInputLines >> szIFlag;
00444 
00445       // error checking
00446       if( (strcmp (szVolAlias.c_str(), "") == 0) ||
00447           (strcmp (szVolId.c_str(), "") == 0))
00448         IOErrorHandler(EPIN);
00449       if( nInputLines < 0 )
00450         IOErrorHandler(ENEGATIVE);
00451 
00452       m_Pincell(i).SetLineOne (szVolId, szVolAlias, nInputLines);
00453       if(szIFlag == "intersect"){
00454         m_Pincell(i).SetIntersectFlag(1);
00455       }
00456       else{
00457         m_Pincell(i).SetIntersectFlag(0);
00458       }
00459       for(int l=1; l<=nInputLines; l++){
00460         if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
00461                                  MAXCHARS, szComment))
00462           IOErrorHandler (INVALIDINPUT);
00463         if (szInputString.substr(0,5) == "pitch"){
00464 
00465           std::istringstream szFormatString (szInputString);
00466           std::cout << "\ngetting pitch data";
00467           szFormatString >> card >> dFlatF >> dLZ;
00468           if( dFlatF < 0 || dLZ < 0  || szFormatString.fail())
00469             IOErrorHandler(ENEGATIVE);
00470           m_Pincell(i).SetPitch (dFlatF, dLZ);
00471         }
00472         if (szInputString.substr(0,9) == "materials"){
00473 
00474           std::istringstream szFormatString (szInputString);
00475           szFormatString >> card >> nMaterials;
00476           if(szFormatString.fail())
00477             IOErrorHandler(INVALIDINPUT);
00478           //setting local arrays
00479           szVMatName.SetSize(nMaterials);
00480           szVMatAlias.SetSize(nMaterials);
00481 
00482           //set class variable sizes
00483           m_Pincell(i).SetMatArray(nMaterials);
00484           std::cout << "\ngetting material data";
00485           for(int j=1; j<= nMaterials; j++){
00486             szFormatString >> szVMatName(j) >> szVMatAlias(j);
00487             if(szFormatString.fail())
00488               IOErrorHandler(INVALIDINPUT);
00489           }
00490           m_Pincell(i).SetMat(szVMatName, szVMatAlias);
00491         }
00492         if (szInputString.substr(0,8) == "cylinder"){
00493 
00494           ++nCyl;
00495           std::cout << "\ngetting cylinder data";
00496           std::istringstream szFormatString (szInputString);
00497           szFormatString >> card >> nRadii >> dVCoor(1) >> dVCoor(2);
00498           if(szFormatString.fail())
00499             IOErrorHandler(INVALIDINPUT);
00500           m_Pincell(i).SetCylSizes(nCyl, nRadii);
00501           m_Pincell(i).SetCylPos(nCyl, dVCoor);
00502 
00503           //set local array
00504           dVCylRadii.SetSize(nRadii);
00505           szVCylMat.SetSize(nRadii);
00506           dVCylZPos.SetSize(2);
00507           //
00508           m_Pincell(i).SetCylSizes(nCyl, nRadii);
00509 
00510           // reading ZCoords - max and min 2 always
00511           for(int k=1; k<=2; k++)
00512             szFormatString >> dVCylZPos(k);
00513           m_Pincell(i).SetCylZPos(nCyl, dVCylZPos);
00514 
00515           // reading Radii
00516           for(int l=1; l<= nRadii; l++){
00517             szFormatString >> dVCylRadii(l);
00518             if( dVCylRadii(l) < 0 || szFormatString.fail())
00519               IOErrorHandler(ENEGATIVE);
00520           }
00521           m_Pincell(i).SetCylRadii(nCyl, dVCylRadii);
00522 
00523           // reading Material alias
00524           for(int m=1; m<= nRadii; m++){
00525             szFormatString >> szVCylMat(m);
00526             if(strcmp (szVCylMat(m).c_str(), "") == 0 || szFormatString.fail())
00527               IOErrorHandler(EALIAS);
00528           }
00529 
00530           m_Pincell(i).SetCylMat(nCyl, szVCylMat);
00531         }
00532         if (szInputString.substr(0,12) == "cellmaterial"){
00533 
00534           std::cout << "\ngetting cell material data";
00535           std::istringstream szFormatString (szInputString);
00536           szFormatString >> card;
00537 
00538           //set local arrays
00539           m_Pincell(i).GetCellMatSize(nCellMat); // since size of cell material already set equal to number of cylinders
00540           dZVStart.SetSize(nCellMat);
00541           dZVEnd.SetSize(nCellMat);
00542           szVCellMat.SetSize(nCellMat);
00543 
00544           for(int k=1; k<=nCellMat; k++){
00545             szFormatString >> dZVStart(k)>> dZVEnd(k) >> szVCellMat(k);
00546             if(strcmp (szVCellMat(k).c_str(), "") == 0 || szFormatString.fail())
00547               IOErrorHandler(EALIAS);
00548           }
00549           m_Pincell(i).SetCellMat(dZVStart, dZVEnd, szVCellMat);
00550         }
00551       }
00552     }// if hexagonal end
00553 
00554   }
00555 
00556 
00557   void AssyGen::ReadAndCreate()
00558   //---------------------------------------------------------------------------
00559   //Function: reads the input file and creates assembly
00560   //Input:    none
00561   //Output:   none
00562   //---------------------------------------------------------------------------
00563   {
00564     //Rewind the input file
00565     m_FileInput.clear (std::ios_base::goodbit);
00566     m_FileInput.seekg (0L, std::ios::beg);
00567     m_nLineNumber = 0;
00568     CParser Parse;
00569     std::string card;
00570 
00571     // start reading the input file break when encounter end
00572     for(;;){
00573       if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
00574                                MAXCHARS, szComment))
00575         IOErrorHandler (INVALIDINPUT);
00576       if (szInputString.substr(0,12) == "geometrytype"){
00577         std::istringstream szFormatString (szInputString);
00578         szFormatString >> card >> m_szGeomType;
00579         if( ((strcmp (m_szGeomType.c_str(), "hexagonal") != 0) &&
00580              (strcmp (m_szGeomType.c_str(), "rectangular") != 0)) || szFormatString.fail())
00581           IOErrorHandler(EGEOMTYPE);
00582 
00583         // set the number of sides in the geometry
00584         if(m_szGeomType == "hexagonal")
00585           m_nSides = 6;
00586         else  if(m_szGeomType == "rectangular")
00587           m_nSides = 4;
00588       }
00589 
00590 
00591       if (szInputString.substr(0,8) == "geometry"){
00592         std::string outfile;
00593         std::istringstream szFormatString (szInputString);
00594         szFormatString >> card >> outfile;
00595         if(strcmp (outfile.c_str(), "surface") == 0 || szFormatString.fail())
00596           m_nPlanar=1;
00597       }
00598       if ((szInputString.substr(0,9) == "materials") && (szInputString.substr(0,19) != "materialset_startid")){
00599 
00600         std::istringstream szFormatString (szInputString);
00601         szFormatString >> card >> m_nAssemblyMat;
00602         if(szFormatString.fail())
00603           IOErrorHandler(INVALIDINPUT);
00604         m_szAssmMat.SetSize(m_nAssemblyMat); m_szAssmMatAlias.SetSize(m_nAssemblyMat);
00605         for (int j=1; j<=m_nAssemblyMat; j++){
00606           szFormatString >> m_szAssmMat(j) >> m_szAssmMatAlias(j);
00607           if( (strcmp (m_szAssmMat(j).c_str(), "") == 0) ||
00608               (strcmp (m_szAssmMatAlias(j).c_str(), "") == 0)){
00609             IOErrorHandler(EMAT);
00610           }
00611           // checking if & inserted at the end of the material by mistake
00612           if (j == m_nAssemblyMat){
00613             std::string dummy = "";
00614             szFormatString >> dummy;
00615             if (strcmp (dummy.c_str(), "") != 0)
00616               IOErrorHandler(EMAT);
00617           }
00618         }
00619       }
00620       if( (szInputString.substr(0,10) == "dimensions") ||
00621           (szInputString.substr(0,4) == "duct") ){
00622 
00623         ++m_nDuctNum;
00624         std::cout << "getting assembly dimensions " << m_nDuctNum << std::endl;
00625 
00626         if(m_szGeomType =="hexagonal"){
00627           std::istringstream szFormatString (szInputString);
00628 
00629           if(m_nDuctNum == 1){
00630             m_dMXYAssm.SetSize(m_nDuct, 2); m_dMZAssm.SetSize(m_nDuct, 2);
00631           }
00632           szFormatString >> card >> m_nDimensions
00633                          >> m_dMXYAssm(m_nDuctNum, 1) >> m_dMXYAssm(m_nDuctNum, 2)
00634                          >> m_dMZAssm(m_nDuctNum, 1) >> m_dMZAssm(m_nDuctNum, 2);
00635           if(m_nDuctNum == 1){
00636             m_dMAssmPitch.SetSize(m_nDuct, m_nDimensions); m_szMMAlias.SetSize(m_nDuct, m_nDimensions);
00637 
00638             assms.resize(m_nDimensions*m_nDuct); // setup while reading the problem size
00639           }
00640 
00641           for (int i=1; i<=m_nDimensions; i++){
00642             szFormatString >> m_dMAssmPitch(m_nDuctNum, i);
00643             if( m_dMAssmPitch(m_nDuctNum, i) < 0 )
00644               IOErrorHandler(ENEGATIVE);
00645           }
00646 
00647           for (int i=1; i<=m_nDimensions; i++){
00648             szFormatString >> m_szMMAlias(m_nDuctNum, i);
00649             if(strcmp (m_szMMAlias(m_nDuctNum, i).c_str(), "") == 0)
00650               IOErrorHandler(EALIAS);
00651           }
00652         }
00653         if(m_szGeomType =="rectangular"){
00654           std::istringstream szFormatString (szInputString);
00655           if(m_nDuctNum == 1){
00656             m_dMXYAssm.SetSize(m_nDuct, 2);
00657             m_dMZAssm.SetSize(m_nDuct, 2);
00658           }
00659           szFormatString >> card >> m_nDimensions
00660                          >> m_dMXYAssm(m_nDuctNum, 1) >> m_dMXYAssm(m_nDuctNum, 2)
00661                          >> m_dMZAssm(m_nDuctNum, 1) >> m_dMZAssm(m_nDuctNum, 2);
00662           if (szFormatString.fail())
00663             IOErrorHandler(INVALIDINPUT);
00664           if(m_nDuctNum == 1){
00665             m_dMAssmPitchX.SetSize(m_nDuct, m_nDimensions);
00666             m_dMAssmPitchY.SetSize(m_nDuct, m_nDimensions);
00667             m_szMMAlias.SetSize(m_nDuct, m_nDimensions);
00668             assms.resize(m_nDimensions*m_nDuct);
00669           }
00670           for (int i=1; i<=m_nDimensions; i++){
00671             szFormatString >> m_dMAssmPitchX(m_nDuctNum, i) >> m_dMAssmPitchY(m_nDuctNum, i);
00672             if( m_dMAssmPitchX(m_nDuctNum, i) < 0 || m_dMAssmPitchY(m_nDuctNum, i) < 0 || szFormatString.fail())
00673               IOErrorHandler(ENEGATIVE);
00674           }
00675 
00676           for (int i=1; i<=m_nDimensions; i++){
00677             szFormatString >> m_szMMAlias(m_nDuctNum, i);
00678             if(strcmp (m_szMMAlias(m_nDuctNum, i).c_str(), "") == 0 || szFormatString.fail())
00679               IOErrorHandler(EALIAS);
00680           }
00681         }
00682       }
00683       if (szInputString.substr(0,8) == "pincells"){
00684         std::istringstream szFormatString (szInputString);
00685 
00686         szFormatString >> card >> m_nPincells >> m_dPitch;
00687         if(m_nPincells < 0)
00688           IOErrorHandler(ENEGATIVE);
00689 
00690         // this is an option if a user wants to specify pitch here
00691         double dTotalHeight = 0.0;
00692 
00693         //get the number of cylinder in each pincell
00694         int nTemp = 1;
00695         if(m_nDimensions > 0){
00696           dTotalHeight = m_dMZAssm(nTemp, 2)-m_dMZAssm(nTemp, 1);
00697         }
00698         else{
00699           dTotalHeight = 0; // nothing specified only pincells in the model
00700         }
00701 
00702         // loop thro' the pincells and read/store pincell data
00703         for (int i=1; i<=m_nPincells; i++){
00704 
00705           // set pitch if specified in pincell card
00706           if(m_dPitch > 0.0)
00707             m_Pincell(i).SetPitch(m_dPitch, dTotalHeight);
00708 
00709           ReadPinCellData( i);
00710           std::cout << "\nread pincell " << i << std::endl;
00711         }
00712       }
00713       if (szInputString.substr(0,8) == "assembly"){
00714         if(m_szGeomType =="hexagonal"){
00715           Create_HexAssm( szInputString);
00716         }
00717         if(m_szGeomType =="rectangular"){
00718           Create_CartAssm(szInputString);
00719         }
00720         if (m_nJouFlag == 0){
00721           CreateOuterCovering();
00722 
00723           // subtract pins before save
00724           Subtract_Pins();
00725           if(m_nPlanar ==1){
00726             Create2DSurf();
00727           }
00728         }
00729       }
00730 
00731       // section the assembly as described in section card
00732       if (szInputString.substr(0,7) == "section" && m_nJouFlag == 0){
00733         std::cout << "Sectioning geometry .." << std::endl;
00734         char cDir;
00735         double dOffset;
00736         std::string szReverse = "";
00737         std::istringstream szFormatString (szInputString);
00738         szFormatString >> card >> cDir >> dOffset >> szReverse;
00739         Section_Assm( cDir, dOffset, szReverse);
00740         std::cout <<"--------------------------------------------------"<<std::endl;
00741 
00742       }
00743       if (szInputString.substr(0,4) == "move" && m_nJouFlag == 0){
00744         std::cout << "Moving geometry .." << std::endl;
00745         double dX, dY, dZ;
00746         std::istringstream szFormatString (szInputString);
00747         szFormatString >> card >> dX >> dY >> dZ;
00748         if(szFormatString.fail())
00749           IOErrorHandler(INVALIDINPUT);
00750         Move_Assm( dX, dY, dZ);
00751         std::cout <<"--------------------------------------------------"<<std::endl;
00752 
00753       }
00754       // ceter the assembly
00755       if (szInputString.substr(0,6) == "center" && m_nJouFlag == 0){
00756 
00757         char rDir = 'q';
00758         std::istringstream szFormatString (szInputString);
00759         szFormatString >> card >> rDir;
00760         if (rDir != 'q')
00761           std::cout << "Positioning assembly to "<< rDir << " center" << std::endl;
00762         else
00763           std::cout << "Positioning assembly to xy center" << std::endl;
00764         Center_Assm( rDir);
00765         std::cout <<"--------------------------------------------------"<<std::endl;
00766       }
00767       // rotate the assembly if rotate card is specified
00768       if (szInputString.substr(0,6) == "rotate" && m_nJouFlag == 0){
00769         char cDir;
00770         double dAngle;
00771         std::cout << "Rotating geometry .." << std::endl;
00772         std::istringstream szFormatString (szInputString);
00773         szFormatString >> card >> cDir >> dAngle;
00774         if(szFormatString.fail())
00775           IOErrorHandler(INVALIDINPUT);
00776         Rotate_Assm( cDir, dAngle);
00777         std::cout <<"--------------------------------------------------"<<std::endl;
00778 
00779       }
00780       // 'yes' or 'no' for creating sidesets
00781       if (szInputString.substr(0,13) == "createsideset"){
00782         std::istringstream szFormatString (szInputString);
00783         szFormatString >> card >> m_szSideset;
00784         std::cout <<"--------------------------------------------------"<<std::endl;
00785       }
00786       // specify a merge tolerance value for cubit journal file
00787       if (szInputString.substr(0,14) == "mergetolerance"){
00788         std::istringstream szFormatString (szInputString);
00789         szFormatString >> card >> m_dMergeTol;
00790         std::cout <<"--------------------------------------------------"<<std::endl;
00791       }  // Handle mesh size inputs
00792       if (szInputString.substr(0,14) == "radialmeshsize"){
00793         std::istringstream szFormatString (szInputString);
00794         szFormatString >> card >> m_dRadialSize;
00795         if(m_dRadialSize < 0 || szFormatString.fail())
00796           IOErrorHandler(ENEGATIVE);
00797         std::cout <<"--------------------------------------------------"<<std::endl;
00798 
00799       }
00800       // Handle mesh size inputs
00801       if (szInputString.substr(0,11) == "tetmeshsize"){
00802         std::istringstream szFormatString (szInputString);
00803         szFormatString >> card >> m_dTetMeshSize;
00804         if(m_dTetMeshSize < 0 || szFormatString.fail())
00805           IOErrorHandler(ENEGATIVE);
00806         std::cout <<"--------------------------------------------------"<<std::endl;
00807 
00808       }
00809       // Handle mesh size inputs
00810       if (szInputString.substr(0,13) == "axialmeshsize"){
00811         std::istringstream szFormatString (szInputString);
00812         szFormatString >> card >> m_dAxialSize;
00813         if(m_dAxialSize < 0 || szFormatString.fail())
00814           IOErrorHandler(ENEGATIVE);
00815         std::cout <<"--------------------------------------------------"<<std::endl;
00816 
00817       }
00818       // Handle mesh size inputs
00819       if (szInputString.substr(0,18) == "neumannset_startid"){
00820         std::istringstream szFormatString (szInputString);
00821         szFormatString >> card >> m_nNeumannSetId;
00822         if(m_nNeumannSetId < 0 || szFormatString.fail())
00823           IOErrorHandler(ENEGATIVE);
00824         std::cout <<"--------------------------------------------------"<<std::endl;
00825 
00826       }
00827       // Handle mesh size inputs
00828       if (szInputString.substr(0,19) == "materialset_startid"){
00829         std::istringstream szFormatString (szInputString);
00830         szFormatString >> card >> m_nMaterialSetId;
00831         if(m_nMaterialSetId < 0 || szFormatString.fail())
00832           IOErrorHandler(ENEGATIVE);
00833         std::cout <<"--------------------------------------------------"<<std::endl;
00834 
00835       }
00836       if (szInputString.substr(0,3) == "end"){
00837 
00838 
00839         if ( m_nJouFlag == 0){
00840           // impring merge before saving
00841           // Imprint_Merge();
00842 
00843           // save .sat file
00844           IBERRCHK(igeomImpl->save(m_szGeomFile.c_str()), *igeomImpl);
00845           std::cout << "Normal Termination.\n"<< "Geometry file: " << m_szGeomFile << " saved." << std::endl;
00846         }
00847         break;
00848       }
00849     }
00850 
00851   }
00852 
00853   void AssyGen::CreateCubitJournal()
00854   //---------------------------------------------------------------------------
00855   //Function: Create Cubit Journal File for generating mesh
00856   //Input:    none
00857   //Output:   none
00858   //---------------------------------------------------------------------------
00859   {
00860     // variables
00861     int nColor;
00862     std::string color[21] = {" ", "thistle", "grey", "deepskyblue", "red", "purple",  "green",
00863                              "yellow", "royalblue", "magenta", "cyan", "lightsalmon", "springgreen",
00864                              "gold", "orange", "brown", "pink", "khaki", "black", "aquamurine", "mediumslateblue"};
00865 
00866     // if creating only journal file load the geometry file
00867     if(m_nJouFlag == 1){
00868       IBERRCHK(igeomImpl->load(m_szGeomFile.c_str()), *igeomImpl);
00869     }
00870 
00871     // get the max and min coordinates of the geometry
00872     double x1, y1, z1, x2, y2, z2;
00873     IBERRCHK(igeomImpl->getBoundBox(x1, y1, z1, x2, y2, z2), *igeomImpl);
00874 
00875     int nSideset=m_nNeumannSetId;
00876     std::string szGrp, szBlock, szSurfTop, szSurfBot, szSize, szSurfSide;
00877     double dHeight = 0.0, dMid = 0.0;
00878     int nTemp = 1;
00879     if(m_nDimensions > 0){
00880       dHeight= fabs(z2 - z1);
00881       dMid = z2 - dHeight/2.0;
00882     }
00883 
00884     // writing to template.jou
00885     m_SchemesFile << "## This file is created by rgg program in MeshKit ##\n";
00886     m_SchemesFile << "##Schemes " << std::endl  ;
00887     m_SchemesFile << "#{CIRCLE =\"circle interval 1 fraction 0.8\"}" << std::endl;
00888     m_SchemesFile << "#{HOLE = \"hole rad_interval 2 bias 0.0\"}" << std::endl;
00889     m_SchemesFile << "#{PAVE = \"pave\"}" << std::endl;
00890     m_SchemesFile << "#{MAP = \"map\"}" << std::endl;
00891     m_SchemesFile << "#{SWEEP = \"sweep\"}" << std::endl;
00892     m_SchemesFile << "#{TET = \"tetmesh\"}" << std::endl;
00893     m_SchemesFile << "## Dimensions" << std::endl;
00894     if(m_szGeomType == "hexagonal"){
00895       if(m_nDimensions > 0){
00896         m_SchemesFile << "#{PITCH =" << m_dMAssmPitch(nTemp, m_nDimensions) << "}" << std::endl;
00897       }
00898     }
00899     else if(m_szGeomType == "rectangular"){
00900       if(m_nDimensions > 0){
00901         m_SchemesFile << "#{PITCHX =" << m_dMAssmPitchX(nTemp, m_nDimensions)<< "}" << std::endl;
00902         m_SchemesFile << "#{PITCHY =" << m_dMAssmPitchY(nTemp, m_nDimensions) << "}" << std::endl;
00903       }
00904     }
00905     if( m_nPlanar ==0){
00906       m_SchemesFile << "#{Z_HEIGHT = " << dHeight << "}" << std::endl;
00907       m_SchemesFile << "#{Z_MID = " << dMid << "}" << std::endl;
00908 
00909     }
00910     m_SchemesFile << "##Set Mesh Sizes" << std::endl;
00911 
00912     if (m_szMeshType == "hex"){
00913       // volume only
00914       if(m_nPlanar == 0 ){
00915         if (-1.0 == m_dAxialSize){
00916           m_SchemesFile << "#{AXIAL_MESH_SIZE = 0.1*Z_HEIGHT}" << std::endl;
00917         }
00918         else {
00919           m_SchemesFile << "#{AXIAL_MESH_SIZE = " << m_dAxialSize << "}" << std::endl;
00920         }
00921 
00922         // create templates for specifying block z intervals
00923         if (m_nDuct > 1){
00924           m_SchemesFile << "## Set interval along Z direction ## " << std::endl;
00925 
00926           for( int p=1; p<= m_nDuct; p++){
00927             m_SchemesFile << "#{BLOCK" << p << "_Z_INTERVAL = AXIAL_MESH_SIZE}" << std::endl;
00928           }
00929           m_SchemesFile << "##" << std::endl;
00930         }
00931       }
00932       if (-1.0 == m_dRadialSize) {
00933         if (m_szGeomType == "hexagonal")
00934           m_SchemesFile << "#{RADIAL_MESH_SIZE = 0.1*PITCH}" << std::endl;
00935         else
00936           m_SchemesFile << "#{RADIAL_MESH_SIZE = 0.02*0.5*(PITCHX+PITCHY)}" << std::endl;
00937       }
00938       else
00939         m_SchemesFile << "#{RADIAL_MESH_SIZE = " << m_dRadialSize << "}" << std::endl;
00940     }
00941     else if (m_szMeshType == "tet"){
00942       if (-1.0 == m_dTetMeshSize) {
00943         if (m_szGeomType == "hexagonal")
00944           m_SchemesFile << "#{TET_MESH_SIZE = 0.1*PITCH}" << std::endl;
00945         else
00946           m_SchemesFile << "#{TET_MESH_SIZE = 0.02*0.5*(PITCHX+PITCHY)}" << std::endl;
00947       }
00948       else {
00949         m_SchemesFile << "#{TET_MESH_SIZE = " << m_dTetMeshSize  << "}" << std::endl;
00950       }
00951     }
00952     // writing schemes .jou file ends, now write the main journal file.
00953 
00954 
00955     // stuff common to both surface and volume
00956     m_FileOutput << "## This file is created by rgg program in MeshKit ##\n";
00957     m_FileOutput << "#User needs to specify mesh interval and schemes in this file\n#" << std::endl;
00958     m_FileOutput << "{include(\"" << m_szSchFile << "\")}" <<std::endl;
00959     m_FileOutput << "#" << std::endl;
00960 
00961     // import the geometry file
00962     m_FileOutput << "# Import geometry file " << std::endl;
00963     m_FileOutput << "import '" << m_szGeomFile <<"'" <<std::endl;
00964     m_FileOutput << "#" << std::endl;
00965 
00966     if(m_szMeshType == "hex"){
00967       // imprint
00968       m_FileOutput << "Merge Tolerance " << m_dMergeTol << std::endl;
00969       m_FileOutput << "#" << std::endl;
00970       m_FileOutput << "#Imprint geometry" << std::endl;
00971       m_FileOutput << "imprint all" << std::endl;
00972       m_FileOutput << "#" << std::endl;
00973 
00974       // merge
00975       m_FileOutput << "#Merge geometry" << std::endl;
00976       m_FileOutput << "merge all" << std::endl;
00977       m_FileOutput << "#" << std::endl;
00978     }
00979 
00980     if(m_szSideset == "yes"){
00981       // top surface sidesets
00982       m_FileOutput << "#Creating top surface sidesets" << std::endl;
00983       for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
00984         ++nSideset;
00985         szSurfTop = m_szAssmMat(p)+"_top";
00986         m_FileOutput << "group 'tmpgrp' equals surface name \""  << szSurfTop  << "\"" << std::endl;
00987         m_FileOutput << "sideset " << nSideset << " surface in tmpgrp" << std::endl;
00988       }
00989       m_FileOutput << "#" << std::endl;
00990     }
00991 
00992     //surface only
00993     if(m_nPlanar ==1){
00994       m_FileOutput << "# Pointing surface normals to 0.0, 0.0, -1.0 or -ve Z or correct STAR-CCM+ cell-face orientation" << std::endl;
00995       m_FileOutput << "surface all normal opposite" << std::endl;
00996       if(m_szSideset == "yes"){
00997         for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
00998           ++nSideset;
00999           szSurfTop = m_szAssmMat(p)+"_top";
01000           m_FileOutput <<  "group 'tmpgrp' equals surface name \""  << szSurfTop  << "\"" << std::endl;
01001           m_FileOutput <<"group 'tmp1' equals curve in surface in tmpgrp" << std::endl;
01002           m_FileOutput << "sideset " << nSideset << " curve in tmp1" << std::endl;
01003         }
01004       }
01005       // group creation dumps. each material surface  has a group
01006       m_FileOutput << "#Creating groups" << std::endl;
01007       for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
01008         szGrp = "g_"+ m_szAssmMat(p);
01009         m_szAssmMat(p);
01010         m_FileOutput << "group \"" << szGrp << "\" add surface name \"" << m_szAssmMat(p) <<"\"" << std::endl;
01011       }
01012       m_FileOutput << "#" << std::endl;
01013 
01014       // block creation dumps
01015       m_FileOutput << "#Creating blocks, Note: you might need to combine some blocks" << std::endl;
01016       for(int p=1; p <= m_szAssmMatAlias.GetSize();p++){
01017         szBlock = "b_"+ m_szAssmMat(p);
01018         szGrp = "g_"+ m_szAssmMat(p);
01019         m_FileOutput << "block " << m_nMaterialSetId + p << " surface in " << szGrp  << std::endl;
01020         m_FileOutput << "block " << m_nMaterialSetId + p << " name \"" << szBlock <<"\""<< std::endl;
01021       }
01022       m_FileOutput << "#" << std::endl;
01023     }
01024 
01025     // volume only
01026     else{
01027       if(m_szSideset == "yes"){
01028         // bottom surface sidesets
01029         m_FileOutput << "#Creating bot and side surface sidesets" << std::endl;
01030 
01031         // rename the skin surfaces, so that they don't appear as sidesets
01032 
01033         for (int p=1; p<=m_nDuct; p++){
01034           for(int q=1;q<=m_nSides; q++){
01035             m_FileOutput << "group 'edge" << (m_nSides*(p-1) + q ) <<"' equals curve with name 'side_edge"
01036                          << (m_nSides*(p-1) + q ) << "@'" << std::endl;
01037 
01038             m_FileOutput << "group 'vt" <<  (m_nSides*(p-1) + q )  <<"' equals vertex with z_max == z_min in curve in edge"
01039                          <<  (m_nSides*(p-1) + q ) << std::endl;
01040 
01041           }
01042         }
01043 
01044 
01045         // creating groups for vertices on the top surface of the duct
01046         for (int p=1; p<=m_nDuct; p++){
01047           for(int q=1;q<=m_nSides; q++){
01048 
01049             if(q != m_nSides){
01050               m_FileOutput << "group 'v" << (m_nSides*(p-1) + q ) <<"' intersect group vt" << (m_nSides*(p-1) + q )
01051                            << " with group vt" <<  (m_nSides*(p-1) + q + 1 )  << std::endl;
01052             }
01053             else {
01054               m_FileOutput << "group 'v" << (m_nSides*(p-1) + q ) <<"' intersect group vt" << (m_nSides*(p-1) + q )
01055                            << " with group vt" <<  (m_nSides*(p-1) + 1 )  << std::endl;
01056             }
01057           }
01058         }
01059         // creating temp surfaces groups
01060         for (int p=1; p<=m_nDuct; p++){
01061           for(int q=1;q<=m_nSides; q++){
01062             m_FileOutput << "group 'st" << (m_nSides*(p-1) + q ) <<"' equals surface with z_max <> z_min in vert in v"
01063                          << (m_nSides*(p-1) + q ) << "'" << std::endl;
01064           }
01065         }
01066 
01067         // creating surface groups for obtaining surfaces
01068         for (int p=1; p<=m_nDuct; p++){
01069           for(int q=1;q<=m_nSides; q++){
01070             if(q != 1){
01071               m_FileOutput << "group 's" << (m_nSides*(p-1) + q ) <<"' intersect group st"  << (m_nSides*(p-1) + q )
01072                            << " with group st" <<  (m_nSides*(p-1) + q - 1 )  << std::endl;
01073             }
01074             else {
01075               m_FileOutput << "group 's" << (m_nSides*(p-1) + q ) <<"' intersect group st" << (m_nSides*(p-1) + q )
01076                            << " with group st" <<  (m_nSides*(p-1) + m_nSides )  << std::endl;
01077             }
01078           }
01079         }
01080 
01081         // renaming the skin side surfaces
01082         for (int p=1; p<=m_nDuct; p++){
01083           for(int q=1;q<=m_nSides; q++){
01084             m_FileOutput << "surface in group s" <<  (m_nSides*(p-1) + q ) << " rename 'side_surface"
01085                          <<  (m_nSides*(p-1) + q ) << "'" << std::endl;
01086 
01087           }
01088         }
01089         if(m_szSideset == "yes"){
01090           // now create top and bot sideset
01091           for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
01092             szSurfTop = m_szAssmMat(p)+"_bot";
01093             szSurfSide = m_szAssmMat(p)+"_side";
01094 
01095 
01096             m_FileOutput << "#" << std::endl;
01097 
01098             ++nSideset;
01099             m_FileOutput << "group 'tmpgrp' equals surface name \""  << szSurfTop  << "\"" << std::endl;
01100             m_FileOutput << "sideset " << nSideset << " surface in tmpgrp" << std::endl;
01101 
01102             ++nSideset;
01103             m_FileOutput << "group 'tmpgrp' equals surface name \""  << szSurfSide  << "\"" << std::endl;
01104             m_FileOutput << "sideset " << nSideset << " surface in tmpgrp" << std::endl;
01105           }
01106           m_FileOutput << "#" << std::endl;
01107         }
01108       }
01109       // group creation dumps. each material surface  has a group
01110       for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
01111         szGrp = "g_"+ m_szAssmMat(p);
01112         m_szAssmMat(p);
01113         m_FileOutput << "group \"" << szGrp << "\" add body name \"" << m_szAssmMat(p) <<"\"" << std::endl;
01114       }
01115       m_FileOutput << "#" << std::endl;
01116 
01117       // block creation dumps
01118       m_FileOutput << "#Creating blocks, Note: you might need to combine some blocks" << std::endl;
01119       for(int p = 1; p <=  m_szAssmMatAlias.GetSize();p++){
01120         szBlock = "b_"+ m_szAssmMat(p);
01121         szGrp = "g_"+ m_szAssmMat(p);
01122         m_FileOutput << "block " <<  m_nMaterialSetId + p << " body in " << szGrp  << std::endl;
01123         m_FileOutput << "block " << m_nMaterialSetId + p << " name \"" << szBlock <<"\""<< std::endl;
01124       }
01125       m_FileOutput << "#" << std::endl;
01126 
01127       if(m_szMeshType == "hex"){
01128 
01129         //now set the sizes
01130         m_FileOutput << "#Set Meshing Scheme and Sizes, use template.jou to specify sizes" << std::endl;
01131 
01132         for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
01133           szGrp = "g_"+ m_szAssmMat(p);
01134           szSize =  m_szAssmMat(p) + "_size";
01135           szSurfBot = m_szAssmMat(p) + "_bot";
01136           szSize =  m_szAssmMat(p) + "_surf_size";
01137           m_FileOutput << "group 'tmpgrp' equals surface name \""  << szSurfBot  << "\"" << std::endl;
01138           m_FileOutput << "surface in tmpgrp  size {"  << szSize <<"}" << std::endl;
01139         }
01140         m_FileOutput << "#" << std::endl;
01141       }
01142     }
01143 
01144     if(m_szMeshType == "hex"){
01145       // some more common stuff meshing top surfaces set the sizes and mesh
01146       m_FileOutput << "#Surfaces mesh, use template.jou to specify sizes" << std::endl;
01147       for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
01148         szSurfTop = m_szAssmMat(p) + "_top";
01149         szGrp = "g_"+ m_szAssmMat(p);
01150         szSize =  m_szAssmMat(p) + "_surf_size";
01151         m_FileOutput << "group 'tmpgrp' equals surface name \""  << szSurfTop  << "\"" << std::endl;
01152         m_FileOutput << "surface in tmpgrp  size {"  << szSize <<"}" << std::endl;
01153         m_FileOutput << "surface in tmpgrp scheme {" << "PAVE" << "}"  << std::endl;
01154         //    m_FileOutput << "mesh surface in " << szGrp << "\n#" << std::endl;
01155 
01156         // dumping these sizes schemes.jou also
01157         m_SchemesFile << "#{"  << szSize <<" = RADIAL_MESH_SIZE}" << std::endl;
01158       }
01159       m_FileOutput << "#" << std::endl;
01160 
01161       // mesh all command after meshing surface
01162       if (m_nDuct <= 1 ){
01163         m_FileOutput << "group 'tmpgrp' add surface name '_top'" << std::endl;
01164         m_FileOutput << "mesh tmpgrp" << std::endl;
01165       }
01166       else {
01167         m_FileOutput << "#Meshing top surface" << std::endl;
01168         m_FileOutput << "mesh surface with z_coord = " << z2 << std::endl;
01169       }
01170 
01171       if(m_nPlanar == 0){ // volumes only
01172         if (m_nDuct == 1){
01173           m_FileOutput << "surface with z_coord > {Z_MID -.1*Z_HEIGHT}" <<
01174             " and z_coord < {Z_MID + .1*Z_HEIGHT} size {AXIAL_MESH_SIZE}" << std::endl ;
01175           m_FileOutput << "mesh vol all" << std::endl;
01176         }
01177         else if (m_nDuct > 1)
01178           m_FileOutput << "### Setting Z intervals on ducts and meshing along Z " << std::endl;
01179         for( int p=m_nDuct; p>= 1; p--){
01180           if(dMid == 0){ // z - centered
01181             m_FileOutput << "surface with z_coord  > " << m_dMZAssm(p, 1) - dHeight/2.0
01182                          << " and z_coord < " << m_dMZAssm(p, 2) - dHeight/2.0 << " interval " << "{BLOCK" << p << "_Z_INTERVAL}" << std::endl;
01183             m_FileOutput << "mesh vol with z_coord  > " << m_dMZAssm(p, 1) - dHeight/2.0
01184                          << " and z_coord < " << m_dMZAssm(p, 2) - dHeight/2.0 << std::endl;
01185           }
01186           else{
01187             m_FileOutput << "surface with z_coord  > " << m_dMZAssm(p, 1)
01188                          << " and z_coord < " << m_dMZAssm(p, 2) << " interval " << "{BLOCK" << p << "_Z_INTERVAL}" << std::endl;
01189             m_FileOutput << "mesh vol with z_coord  > " << m_dMZAssm(p, 1)
01190                          << " and z_coord < " << m_dMZAssm(p, 2) << std::endl;
01191 
01192             m_FileOutput << "##" << std::endl;
01193           }
01194         }
01195       }
01196     }
01197     else if(m_szMeshType == "tet"){
01198       m_FileOutput << "##"<< std::endl;
01199       m_FileOutput << "# groupings for creating vertex groups"<< std::endl;
01200       for (int p=1; p<=m_nDuct; p++){
01201         for(int q=1;q<=m_nSides; q++){
01202           m_FileOutput << "group 'edge" << (m_nSides*(p-1) + q ) <<"' equals curve with name 'side_edge"
01203                        << (m_nSides*(p-1) + q ) << "@'" << std::endl;
01204 
01205           m_FileOutput << "group 'vt" <<  (m_nSides*(p-1) + q )  <<"' equals vertex with z_max == z_min in curve in edge"
01206                        <<  (m_nSides*(p-1) + q ) << std::endl;
01207 
01208         }
01209       }
01210 
01211 
01212       // creating groups for vertices on the top surface of the duct
01213       for (int p=1; p<=m_nDuct; p++){
01214         for(int q=1;q<=m_nSides; q++){
01215 
01216           if(q != m_nSides){
01217             m_FileOutput << "group 'v" << (m_nSides*(p-1) + q ) <<"' intersect group vt" << (m_nSides*(p-1) + q )
01218                          << " with group vt" <<  (m_nSides*(p-1) + q + 1 )  << std::endl;
01219           }
01220           else {
01221             m_FileOutput << "group 'v" << (m_nSides*(p-1) + q ) <<"' intersect group vt" << (m_nSides*(p-1) + q )
01222                          << " with group vt" <<  (m_nSides*(p-1) + 1 )  << std::endl;
01223           }
01224         }
01225       }
01226 
01227       // creating temp surfaces groups
01228       for (int p=1; p<=m_nDuct; p++){
01229         for(int q=1;q<=m_nSides; q++){
01230           m_FileOutput << "group 'st" << (m_nSides*(p-1) + q ) <<"' equals surface with z_max <> z_min in vert in v"
01231                        << (m_nSides*(p-1) + q ) << "'" << std::endl;
01232         }
01233       }
01234 
01235       // creating side curve and surface groups
01236       for (int p=1; p<=m_nDuct; p++){
01237         for(int q=1;q<=m_nSides; q++){
01238 
01239           m_FileOutput << "group 'c" <<  (m_nSides*(p-1) + q )  <<"' equals curve with z_max <> z_min in vert in v"
01240                        <<  (m_nSides*(p-1) + q ) << std::endl;
01241 
01242           if(q != 1){
01243             m_FileOutput << "group 's" << (m_nSides*(p-1) + q ) <<"' intersect group st"  << (m_nSides*(p-1) + q )
01244                          << " with group st" <<  (m_nSides*(p-1) + q - 1 )  << std::endl;
01245           }
01246           else {
01247             m_FileOutput << "group 's" << (m_nSides*(p-1) + q ) <<"' intersect group st" << (m_nSides*(p-1) + q )
01248                          << " with group st" <<  (m_nSides*(p-1) + m_nSides )  << std::endl;
01249           }
01250         }
01251       }
01252 
01253       // renaming the side surfaces for getting the split surfaces later
01254       for (int p=1; p<=m_nDuct; p++){
01255         for(int q=1;q<=m_nSides; q++){
01256           m_FileOutput << "surface in group s" <<  (m_nSides*(p-1) + q ) << " rename 'side_surface"
01257                        <<  (m_nSides*(p-1) + q ) << "'" << std::endl;
01258         }
01259       }
01260 
01261       // splitting the surfaces
01262       for (int p=1; p<=m_nDuct; p++){
01263         for(int q=1;q<=m_nSides; q++){
01264 
01265           m_FileOutput << "split surface in group s" <<  (m_nSides*(p-1) + q )  <<" direction curve in group c"
01266                        <<  (m_nSides*(p-1) + q ) << std::endl;
01267         }
01268       }
01269 
01270       // get all the split surfaces in individual groups
01271       for (int p=1; p<=m_nDuct; p++){
01272         for(int q=1;q<=m_nSides; q++){
01273 
01274           m_FileOutput << "group 'sname" << (m_nSides*(p-1) + q ) <<  "' equals surface with name 'side_surface"
01275                        <<  (m_nSides*(p-1) + q ) << "'"<< std::endl;
01276           m_FileOutput << "group 'svert" << (m_nSides*(p-1) + q ) <<  "' equals surface in vert in v"
01277                        <<  (m_nSides*(p-1) + q ) << std::endl;
01278           m_FileOutput << "group 'ssplit" << (m_nSides*(p-1) + q ) <<  "' intersect group sname" <<  (m_nSides*(p-1) + q )
01279                        << " with group svert" << (m_nSides*(p-1) + q ) << std::endl;
01280         }
01281       }
01282 
01283       // get all the split surfaces in individual groups
01284       for (int p=1; p<=m_nDuct; p++){
01285         for(int q=1;q<=m_nSides; q++){
01286 
01287           if(q != 1){
01288             m_FileOutput << "group 'ssplit_" << (m_nSides*(p-1) + q ) <<"' intersect group sname"  << (m_nSides*(p-1) + q )
01289                          << " with group svert" <<  (m_nSides*(p-1) + q - 1 )  << std::endl;
01290           }
01291           else {
01292             m_FileOutput << "group 'ssplit_" <<  (m_nSides*(p-1) + q ) <<"' intersect group sname" << (m_nSides*(p-1) + q )
01293                          << " with group svert" <<  (m_nSides*(p-1) + m_nSides )  << std::endl;
01294           }
01295         }
01296       }
01297       // imprint
01298       m_FileOutput << "Merge Tolerance " << m_dMergeTol << std::endl;
01299       m_FileOutput << "#" << std::endl;
01300       m_FileOutput << "#Imprint geometry" << std::endl;
01301       m_FileOutput << "imprint all" << std::endl;
01302       m_FileOutput << "#" << std::endl;
01303 
01304       // merge
01305       m_FileOutput << "#Merge geometry" << std::endl;
01306       m_FileOutput << "merge all" << std::endl;
01307       m_FileOutput << "#" << std::endl;
01308 
01309       m_FileOutput << "#Set mesh scheme and size" << std::endl;
01310       m_FileOutput << "volume all scheme {TET} size {TET_MESH_SIZE}" << std::endl;
01311 
01312       // mesh one side of each duct, such that one is flipped mesh of the other
01313       for (int p=1; p<=m_nDuct; p++){
01314 
01315         m_FileOutput << "mesh surface in group ssplit" <<  (m_nSides*(p-1) + 1) << std::endl;
01316 
01317         m_FileOutput << "surface in group ssplit_" << (m_nSides*(p-1) + 1) << " scheme copy source surface in group ssplit"
01318                      <<  (m_nSides*(p-1) + 1)
01319                      << " source  mk = new MKCore(); curve in group c" <<  (m_nSides*(p-1) + 1 ) << " target curve in group c" <<  (m_nSides*(p-1) + m_nSides )
01320                      << " source vertex in group v" <<  (m_nSides*(p-1) + 1) << " target vertex in group v"  <<  (m_nSides*(p-1) + m_nSides )
01321                      << " nosmoothing" << std::endl;
01322 
01323         m_FileOutput << "mesh surface in group  ssplit_" << (m_nSides*(p-1) + 1) << std::endl;
01324       }
01325 
01326       // setting the copy mesh commands on the above pair of split surfaces to have all surfaces symmetrical
01327       for (int p=1; p<=m_nDuct; p++){
01328         for(int q=1;q<=m_nSides; q++){
01329           if(q != m_nSides){
01330             m_FileOutput << "copy mesh surface in ssplit" <<  (m_nSides*(p-1) + 1)
01331                          << " onto surface in ssplit" << (m_nSides*(p-1) + q + 1 )
01332                          << " source vertex in group v" << (m_nSides*(p-1) + 1)
01333                          << " target vertex in group v" <<  (m_nSides*(p-1) + q + 1) << " nosmoothing" <<  std::endl;
01334 
01335             m_FileOutput << "copy mesh surface in ssplit_" << (m_nSides*(p-1) + 1 )
01336                          << " onto surface in ssplit_" << (m_nSides*(p-1) + q +1 )
01337                          << " source vertex in group v" << (m_nSides*p)
01338                          << " target vertex in group v" <<  (m_nSides*(p-1) + q) << " nosmoothing" << std::endl;
01339 
01340           }
01341           else{
01342             // do nothing
01343           }
01344         }
01345       }
01346 
01347       m_FileOutput << "#  mk = new MKCore(); Mesh all volumes now" << std::endl;
01348       m_FileOutput << "mesh vol all" << std::endl;
01349     }
01350     // color now
01351     m_FileOutput << "#Set color for different parts" << std::endl;
01352     if(m_nPlanar == 0){ // volumes only
01353       for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
01354         szGrp = "g_"+ m_szAssmMat(p);
01355         if(p>20)
01356           nColor = 1;
01357         else
01358           nColor = p;
01359         m_FileOutput << "color body in " << szGrp << " " << color[nColor] << std::endl;
01360       }
01361     }
01362     else{ //surfaces
01363       // color now
01364       for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
01365         szGrp = "g_"+ m_szAssmMat(p);
01366         if(p>20)
01367           nColor = 1;
01368         else
01369           nColor = p;
01370         m_FileOutput << "color surface in " << szGrp << " " << color[nColor] << std::endl;
01371       }
01372     }
01373 
01374 
01375     // save as .cub file dump
01376     m_FileOutput << "#\n#Save file" << std::endl;
01377     std::string szSave = m_szFile + ".cub";
01378     m_FileOutput << "save as '"<< szSave <<"'" << " overwrite"<<std::endl;
01379 
01380 
01381     std::cout << "Schemes file created: " << m_szSchFile << std::endl;
01382     std::cout << "Cubit journal file created: " << m_szJouFile << std::endl;
01383 
01384   }
01385 
01386   void AssyGen:: ComputePinCentroid( int nTempPin, CMatrix<std::string> MAssembly,
01387                                      int m, int n, double &dX, double &dY, double &dZ)
01388   // ---------------------------------------------------------------------------
01389   // Function: computes the centroid in the whole assembly of rectangular or hexagonal pincell
01390   // Input:    number and location of the pincell
01391   // Output:   coordinates of pin in assembly
01392   // ---------------------------------------------------------------------------
01393   {
01394     int nTempPin1 = 0, nTempPin2 = 0, nInputLines;
01395     std::string szVolId, szVolAlias;
01396     if(m_szGeomType == "hexagonal"){
01397       double dP, dZ;
01398       m_Pincell(nTempPin).GetPitch(dP, dZ);
01399 
01400       if (m < m_nPin){
01401         dX = (m_nPin - n + 1)*dP/2.0 + n*dP/2.0 + (n-1)*dP - (m-1)*dP/2.0;
01402         dY = (m-1)*(0.5*dP/sin(pi/3.0) + 0.5*dP*sin(pi/6.0)/sin(pi/3.0));
01403       }
01404       else{
01405         dX = (m_nPin - n + 1)*dP/2.0 + n*dP/2.0 + (n-1)*dP - (2*m_nPin - m -1)*dP/2.0;
01406         dY = (m-1)*(0.5*dP/sin(pi/3.0) + 0.5*dP*sin(pi/6.0)/sin(pi/3.0));
01407       }
01408     }
01409     if(m_szGeomType == "rectangular"){
01410       double dPX, dPY, dPZ, dPX1, dPY1, dPZ1, dPX2, dPY2, dPZ2;
01411       if((m_Assembly(m,n)=="x")||(m_Assembly(m,n)=="xx"))
01412         m_Pincell(1).GetPitch(dPX, dPY, dPZ);
01413       else
01414         m_Pincell(nTempPin).GetPitch(dPX, dPY, dPZ);
01415 
01416       if (n==1){
01417         dX = 0;
01418         if(m==1)
01419           dY = 0;
01420       }
01421       else{
01422         dX+= dPX/2.0;
01423         // find the previous pincell type
01424         // check if it's dummy
01425         if((m_Assembly(m,n-1)=="x")||(m_Assembly(m,n-1)=="xx")){
01426           dX+=dPX/2.0;
01427         }
01428         else{
01429           for(int b=1; b<=m_nPincells; b++){
01430             m_Pincell(b).GetLineOne(szVolId, szVolAlias, nInputLines);
01431             if(m_Assembly(m,n-1) == szVolAlias)
01432               nTempPin1 = b;
01433           }
01434           m_Pincell(nTempPin1).GetPitch(dPX1, dPY1, dPZ1);
01435           // now add half of X pitch to the previous cells pitch
01436           dX+= dPX1/2.0;
01437         }
01438       }
01439       if (m > 1 && n==1){
01440         dY+= dPY/2.0;
01441         // check if it's dummy
01442         if((m_Assembly(m-1,n)=="x")||(m_Assembly(m-1,n)=="xx")){
01443           dY+= dPY/2.0;
01444         }
01445         else{
01446           for(int c=1; c<=m_nPincells; c++){
01447             m_Pincell(c).GetLineOne(szVolId, szVolAlias, nInputLines);
01448             if(m_Assembly(m-1,n) == szVolAlias)
01449               nTempPin2 = c;
01450           }
01451           m_Pincell(nTempPin2).GetPitch(dPX2, dPY2, dPZ2);
01452           dY+= dPY2/2.0;
01453         }
01454       }
01455       dZ = 0.0; // moving in XY plane only
01456     }//if rectangular ends
01457 
01458   }
01459 
01460   void AssyGen::IOErrorHandler (ErrorStates ECode) const
01461   // ---------------------------------------------------------------------------
01462   // Function: displays error messages related to input data
01463   // Input:    error code
01464   // Output:   none
01465   // ---------------------------------------------------------------------------
01466   {
01467     std::cerr << '\n';
01468 
01469     if (ECode == PINCELLS) // invalid number of pincells
01470       std::cerr << "Number of pincells must be >= 0.";
01471     else if (ECode == INVALIDINPUT) // invalid input
01472       std::cerr << "Invalid input.";
01473     else if (ECode == EMAT) // invalid input
01474       std::cerr << "Invalid Material Data.";
01475     else if (ECode == EGEOMTYPE) // invalid input
01476       std::cerr << "Invalid GeomType Data.";
01477     else if (ECode == EGEOMENGINE) // invalid input
01478       std::cerr << "Invalid Geometry Engine.";
01479     else if (ECode == EALIAS) // invalid input
01480       std::cerr << "Error Reading Aliases.";
01481     else if (ECode == ENEGATIVE) // invalid input
01482       std::cerr << "Unexpected negative value.";
01483     else if (ECode == EPIN) // invalid input
01484       std::cerr << "Invalid pinCell specs.";
01485     else
01486       std::cerr << "Unknown error ...?";
01487 
01488     std::cerr << '\n' << "Error in input file line : " << m_nLineNumber;
01489     std::cerr << std::endl;
01490     exit (1);
01491   }
01492 
01493   void AssyGen:: Name_Faces( const std::string sMatName, const iBase_EntityHandle body,  iBase_TagHandle this_tag )
01494   // ---------------------------------------------------------------------------
01495   // Function: names all the faces in the body
01496   // Input:    none
01497   // Output:   none
01498   // ---------------------------------------------------------------------------
01499   {
01500     // get the surface with max z
01501     double dTol=1.0e-6, dZTemp = 0.0;
01502     int flag = 0, locTemp = 0;
01503     iBase_EntityHandle max_surf = NULL, min_surf = NULL, side_surf =NULL;
01504     //SimpleArray<iBase_EntityHandle> surfs;
01505     std::vector<iBase_EntityHandle> surfs;
01506     int nSide = 0;
01507     std::ostringstream os;
01508     std::string sMatName0=sMatName+"_top";
01509     std::string sMatName1=sMatName+"_bot";
01510     std::string sSideName;
01511     IBERRCHK(igeomImpl->getEntAdj(body, iBase_FACE, surfs), *igeomImpl);
01512 
01513     //SimpleArray<double> max_corn, min_corn;
01514     //std::vector <double> max_corn, min_corn;
01515 
01516     double *max_corn = new double [3*surfs.size()];
01517     double *min_corn = new double [3*surfs.size()];
01518     IBERRCHK(igeomImpl->getArrBoundBox(&surfs[0], (int) surfs.size(),iBase_INTERLEAVED, &min_corn[0], &max_corn[0]), *igeomImpl);
01519 
01520     for (int i = 0; i < (int) surfs.size(); ++i){
01521       // first find the max z-coordinate
01522       if( (fabs(min_corn[3*i+2]-max_corn[3*i+2])) < dTol ) {
01523         if(flag == 0){
01524           dZTemp = min_corn[3*i+2];
01525           locTemp = i;
01526           flag = 1;
01527         }
01528         else if(dZTemp > min_corn[3*i+2]){
01529           // we have a bot surface
01530           min_surf = surfs[i];
01531           // the top surface is dZTemp
01532           max_surf = surfs[locTemp];
01533         }
01534         else{
01535           //we have a top surface
01536           min_surf = surfs[locTemp];
01537           // the top surface is dZTemp
01538           max_surf = surfs[i];
01539         }
01540       }
01541       // see if max or min set name
01542       if(max_surf !=0){
01543         IBERRCHK(igeomImpl->setData(max_surf, this_tag, sMatName0.c_str()), *igeomImpl);
01544 
01545         std::cout << sMatName0 << ",  ";
01546         max_surf = NULL;
01547 
01548       }
01549       if(min_surf !=0){
01550         IBERRCHK(igeomImpl->setData(min_surf, this_tag, sMatName1.c_str()), *igeomImpl);
01551         std::cout << sMatName1 << ",  ";
01552         min_surf = NULL;
01553 
01554       }
01555     }
01556     for (int i = 0; i < (int) surfs.size(); ++i){
01557       if( (fabs(min_corn[3*i+2]-max_corn[3*i+2])) < dTol ) {
01558         continue; // its a max of min surface
01559       }
01560       else{ // its a side surface
01561         side_surf = surfs[i];
01562       }
01563       //set name for the sidesurf now
01564       if(side_surf !=0){
01565         ++nSide;
01566         sSideName = sMatName + "_side";
01567         os << sSideName << nSide;
01568         sSideName = os.str();
01569         IBERRCHK(igeomImpl->setData(side_surf, this_tag, sMatName1.c_str()), *igeomImpl);
01570         std::cout << sSideName << ",  " ;
01571         sSideName = "";
01572         os.str("");
01573         side_surf = NULL;
01574       }
01575       else {
01576         std::cerr << "Couldn't find surface for naming" << std::endl;
01577       }
01578     }
01579     std::cout <<"\n";
01580 
01581   }
01582 
01583   void AssyGen::Center_Assm ( char &rDir)
01584   // ---------------------------------------------------------------------------
01585   // Function: centers all the entities along x and y axis
01586   // Input:    none
01587   // Output:   none
01588   // ---------------------------------------------------------------------------
01589   {
01590     double xmin, xmax, ymin, ymax, zmin, zmax, xcenter = 0.0, ycenter = 0.0, zcenter = 0.0;
01591     // position the assembly such that origin is at the center before sa
01592     IBERRCHK(igeomImpl->getBoundBox(xmin, ymin,zmin, xmax, ymax, zmax), *igeomImpl);
01593 
01594 
01595     // moving all geom entities to center
01596 
01597 
01598     if( rDir =='x'){
01599       xcenter = (xmin+xmax)/2.0;
01600     }
01601     else if( rDir =='y'){
01602       ycenter = (ymin+ymax)/2.0;
01603     }
01604     else if ( rDir =='z'){
01605       zcenter = (zmin+zmax)/2.0;
01606     }
01607     else{
01608       // assume that it is centered along x and y and not z direction
01609       xcenter = (xmin+xmax)/2.0;
01610       ycenter = (ymin+ymax)/2.0;
01611     }
01612 
01613     //SimpleArray<iBase_EntityHandle> all;
01614     std::vector<iBase_EntityHandle> all;
01615     IBERRCHK(igeomImpl->getEntities(root_set, iBase_REGION, all), *igeomImpl);
01616 
01617 
01618     for(int i=0; i< (int) all.size(); i++){
01619       IBERRCHK(igeomImpl->moveEnt(all[i], -xcenter, -ycenter, -zcenter), *igeomImpl);
01620 
01621     }
01622 
01623   }
01624 
01625   void AssyGen::Section_Assm ( char &cDir, double &dOffset, const std::string szReverse)
01626   // ---------------------------------------------------------------------------
01627   // Function: sections the assembly about the cutting plane
01628   // Input:    none
01629   // Output:   none
01630   // ---------------------------------------------------------------------------
01631   {
01632     double xmin, xmax, ymin, ymax, zmin, zmax, yzplane = 0.0, xzplane = 0.0;
01633     iBase_EntityHandle sec = NULL;
01634     int nReverse = 0;
01635     // check if reverse side is needed
01636     if(szReverse == "reverse"){
01637       nReverse = 1;
01638     }
01639     if( cDir =='x'){
01640       yzplane = 1.0;
01641       xzplane = 0.0;
01642     }
01643     if( cDir =='y'){
01644       yzplane = 0.0;
01645       xzplane = 1.0;
01646     }
01647     //SimpleArray<iBase_EntityHandle> all;
01648     std::vector<iBase_EntityHandle> all;
01649     IBERRCHK(igeomImpl->getEntities(root_set, iBase_REGION, all), *igeomImpl);
01650 
01651     // loop and section/delete entities
01652     for(int i=0; i < (int) all.size(); i++){
01653       //get the bounding box to decide
01654       IBERRCHK(igeomImpl->getEntBoundBox(all[i], xmin, ymin, zmin, xmax, ymax, zmax), *igeomImpl);
01655 
01656       if(xmin > dOffset && yzplane ==1 && nReverse ==1){
01657         IBERRCHK(igeomImpl->deleteEnt(all[i]), *igeomImpl);
01658 
01659         continue;
01660       }
01661       if(ymin > dOffset && xzplane == 1 && nReverse ==1){
01662         IBERRCHK(igeomImpl->deleteEnt(all[i]), *igeomImpl);
01663 
01664         continue;
01665       }
01666       if(xmax < dOffset && yzplane ==1 && nReverse ==0){
01667         IBERRCHK(igeomImpl->deleteEnt(all[i]), *igeomImpl);
01668 
01669         continue;
01670       }
01671       if(ymax < dOffset && xzplane == 1 && nReverse ==0){
01672         IBERRCHK(igeomImpl->deleteEnt(all[i]), *igeomImpl);
01673 
01674         continue;
01675       }
01676       else{
01677         if(xzplane ==1 && ymax >dOffset && ymin < dOffset){
01678           IBERRCHK(igeomImpl->sectionEnt(all[i],yzplane,xzplane,0, dOffset, nReverse, sec), *igeomImpl);
01679 
01680         }
01681         if(yzplane ==1 && xmax >dOffset && xmin < dOffset){
01682           IBERRCHK(igeomImpl->sectionEnt(all[i],yzplane,xzplane,0, dOffset, nReverse, sec), *igeomImpl);
01683 
01684         }
01685       }
01686     }
01687 
01688   }
01689 
01690   void AssyGen::Rotate_Assm ( char &cDir, double &dAngle)
01691   // ---------------------------------------------------------------------------
01692   // Function: rotates the whole assembly
01693   // Input:    none
01694   // Output:   none
01695   // ---------------------------------------------------------------------------
01696   {
01697     double dX = 0.0, dY=0.0, dZ=0.0;
01698     if( cDir =='x'){
01699       dX = 1.0;
01700     }
01701     if( cDir =='y'){
01702       dY = 1.0;
01703     }
01704     if( cDir =='z'){
01705       dZ = 1.0;
01706     }
01707     std::vector<iBase_EntityHandle> all;
01708     IBERRCHK(igeomImpl->getEntities(root_set, iBase_REGION, all), *igeomImpl);
01709 
01710 
01711     // loop and rotate all entities
01712     for(int i=0; i< (int) all.size(); i++){
01713       //get the bounding box to decide
01714       IBERRCHK(igeomImpl->rotateEnt(all[i], dAngle, dX, dY, dZ), *igeomImpl);
01715 
01716     }
01717 
01718   }
01719 
01720   void AssyGen::Move_Assm ( double &dX,double &dY, double &dZ)
01721   // ---------------------------------------------------------------------------
01722   // Function: move's the model by dX, dY and dZ
01723   // Input:    none
01724   // Output:   none
01725   // ---------------------------------------------------------------------------
01726   {
01727     std::vector<iBase_EntityHandle> all;
01728     IBERRCHK(igeomImpl->getEntities(root_set, iBase_REGION, all), *igeomImpl);
01729 
01730 
01731     // loop and rotate all entities
01732     for(int i=0; i< (int) all.size(); i++){
01733       //get the bounding box to decide
01734       IBERRCHK(igeomImpl->moveEnt(all[i], dX, dY, dZ), *igeomImpl);
01735 
01736     }
01737 
01738   }
01739 
01740   void AssyGen::Create_HexAssm( std::string &szInputString)
01741   // ---------------------------------------------------------------------------
01742   // Function: read and create the assembly for hexagonal lattice
01743   // Input:    error code
01744   // Output:   none
01745   // ---------------------------------------------------------------------------
01746   {
01747     CParser Parse;
01748     std::string card, szVolId, szVolAlias;
01749     int nInputLines, nTempPin = 1, t, nIFlag = 0;
01750     double dX = 0.0, dY =0.0, dZ=0.0;
01751     double  dP, dH, dSide, dHeight;
01752     iBase_EntityHandle assm = NULL;
01753     std::cout << "\ngetting Assembly data and creating ..\n"<< std::endl;
01754     std::istringstream szFormatString (szInputString);
01755 
01756     szFormatString >> card >> m_nPin;
01757 
01758     // width of the hexagon is n*n+1/2
01759     int nWidth =2*m_nPin -1;
01760 
01761     // creating a square array of size width
01762     m_Assembly.SetSize(nWidth, nWidth);
01763     if (m_nJouFlag == 1)
01764       return;
01765 
01766     for(int m=1; m<=nWidth; m++){
01767       if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
01768                                MAXCHARS, szComment))
01769         IOErrorHandler (INVALIDINPUT);
01770       if(m>m_nPin)
01771         t = 2*m_nPin - m;
01772       else
01773         t = m;
01774       std::istringstream szFormatString1 (szInputString);
01775 
01776       for(int n=1; n<=(m_nPin + t - 1); n++){
01777         szFormatString1 >> m_Assembly(m,n);
01778         if(szFormatString1.fail())
01779           IOErrorHandler (INVALIDINPUT);
01780         // if dummy pincell skip and continue
01781         if((m_Assembly(m,n)=="x")||(m_Assembly(m,n)=="xx")){
01782           continue;
01783         }
01784         // find that pincell
01785         for(int b=1; b<=m_nPincells; b++){
01786           m_Pincell(b).GetLineOne(szVolId, szVolAlias, nInputLines);
01787           if(m_Assembly(m,n) == szVolAlias)
01788             nTempPin = b;
01789         }
01790 
01791         // now compute the location and create it
01792         ComputePinCentroid( nTempPin, m_Assembly, m, n, dX, dY, dZ);
01793 
01794         // now create the pincell in the location found
01795         std::cout << "\n--------------------------------------------------"<<std::endl;
01796         std::cout << " m = " << m <<" n = " << n << std::endl;
01797         std::cout << "creating pin: " << nTempPin;
01798         std::cout << " at X Y Z " << dX << " " << dY << " " << dZ << std::endl;
01799 
01800         m_Pincell(nTempPin).GetIntersectFlag(nIFlag);
01801         if(nIFlag){
01802           CreatePinCell_Intersect( nTempPin, dX, -dY, dZ);
01803         }
01804         else{
01805           CreatePinCell( nTempPin, dX, -dY, dZ);
01806         }
01807       }
01808     }
01809 
01810     // get all the entities (in pins)defined so far, in an entity set - for subtraction later
01811     IBERRCHK(igeomImpl->getEntities(root_set, iBase_REGION, in_pins), *igeomImpl);
01812 
01813     std::cout << "\n\nCreating surrounding outer hexes .." << std::endl;
01814 
01815     for (int nTemp = 1; nTemp <= m_nDuct; nTemp ++){
01816       if(m_nDimensions >0){
01817 
01818         // create outermost hexes
01819         for(int n=1;n<=m_nDimensions; n++){
01820           dSide = m_dMAssmPitch(nTemp, n)/(sqrt(3));
01821           dHeight = m_dMZAssm(nTemp, 2) - m_dMZAssm(nTemp, 1);
01822 
01823           // creating coverings
01824           IBERRCHK(igeomImpl->createPrism(dHeight, 6, dSide, dSide, assm), *igeomImpl);
01825 
01826 
01827           // rotate the prism to match the pins
01828           IBERRCHK(igeomImpl->rotateEnt(assm, 30, 0, 0, 1), *igeomImpl);
01829 
01830 
01831           m_Pincell(1).GetPitch(dP, dH);
01832           dX = m_nPin*dP;
01833           dY = -(m_nPin-1)*dP*sqrt(3.0)/2.0;
01834           dZ = (m_dMZAssm(nTemp, 2) + m_dMZAssm(nTemp, 1))/2.0;
01835 
01836           // position the prism
01837           IBERRCHK(igeomImpl->moveEnt(assm, dX, dY, dZ), *igeomImpl);
01838 
01839 
01840           // populate the coverings array
01841           assms[(nTemp-1)*m_nDimensions + n -1]=assm;
01842         }
01843       }
01844     }
01845 
01846   }
01847 
01848   void AssyGen::Create_CartAssm( std::string &szInputString)
01849   // ---------------------------------------------------------------------------
01850   // Function: read and create the assembly for rectangular lattice
01851   // Input:    error code
01852   // Output:   none
01853   // ---------------------------------------------------------------------------
01854   {
01855     CParser Parse;
01856     std::string card, szVolId, szVolAlias;
01857     int nInputLines, nTempPin = 1, nIFlag = 0.0;
01858     double dX = 0.0, dY =0.0, dZ=0.0, dMoveX = 0.0, dMoveY = 0.0, dHeight = 0, dPX=0.0, dPY=0.0, dPZ=0.0;
01859     iBase_EntityHandle assm = NULL;
01860 
01861     std::istringstream szFormatString (szInputString);
01862     szFormatString >> card >> m_nPinX >> m_nPinY;
01863     m_Assembly.SetSize(m_nPinX,m_nPinY);
01864 
01865     if (m_nJouFlag == 1)
01866       return;
01867 
01868 
01869     //read the next line to get assembly info &store assembly info
01870     for(int m=1; m<=m_nPinY; m++){
01871       if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
01872                                MAXCHARS, szComment))
01873         IOErrorHandler (INVALIDINPUT);
01874       std::istringstream szFormatString1 (szInputString);
01875 
01876       //store the line read in Assembly array and create / position the pin in the core
01877       for(int n=1; n<=m_nPinX; n++){
01878         szFormatString1 >> m_Assembly(m,n);
01879         if(szFormatString1.fail())
01880           IOErrorHandler (INVALIDINPUT);
01881 
01882 
01883         // loop thro' all pins to get the type of pin
01884         for(int b=1; b<=m_nPincells; b++){
01885           m_Pincell(b).GetLineOne(szVolId, szVolAlias, nInputLines);
01886           if(m_Assembly(m,n) == szVolAlias)
01887             nTempPin = b;
01888         }
01889 
01890         //now compute the location where the pin needs to be placed
01891         ComputePinCentroid( nTempPin, m_Assembly, m, n, dX, dY, dZ);
01892 
01893         // if dummy pincell skip and continue
01894         if((m_Assembly(m,n)=="x")||(m_Assembly(m,n)=="xx")){
01895           m_Pincell(1).GetPitch(dPX, dPY, dPZ);
01896           continue;
01897         }
01898 
01899         // now create the pincell in the location found
01900         std::cout << "\n--------------------------------------------------"<<std::endl;
01901         std::cout << " m = " << m <<" n = " << n << std::endl;
01902         std::cout << "creating pin: " << nTempPin;
01903         std::cout << " at X Y Z " << dX << " " << -dY << " " << dZ << std::endl;
01904 
01905         m_Pincell(nTempPin).GetIntersectFlag(nIFlag);
01906         if(nIFlag){
01907           CreatePinCell_Intersect( nTempPin, dX, -dY, dZ);
01908         }
01909         else{
01910           CreatePinCell( nTempPin, dX, -dY, dZ);
01911         }
01912         // dMoveX and dMoveY are stored for positioning the outer squares later
01913         if(m == m_nPinY && n ==m_nPinX){
01914           dMoveX = dX/2.0;
01915           dMoveY = -dY/2.0;
01916         }
01917       }
01918     }
01919     std::cout << "\n--------------------------------------------------"<<std::endl;
01920 
01921     // get all the entities (in pins)defined so far, in an entity set - for subtraction later
01922     IBERRCHK(igeomImpl->getEntities(root_set, iBase_REGION, in_pins), *igeomImpl);
01923 
01924 
01925 
01926     if(m_nDimensions > 0){
01927 
01928       // create outermost rectangular blocks
01929       std::cout << "\nCreating surrounding outer blocks .." << std::endl;
01930       int nCount = -1;
01931       for(int nTemp = 1; nTemp <= m_nDuct; nTemp ++){
01932         for(int n=1;n<=m_nDimensions; n++){
01933           ++nCount;
01934           dHeight = m_dMZAssm(nTemp, 2) - m_dMZAssm(nTemp, 1);
01935           IBERRCHK(igeomImpl->createBrick(m_dMAssmPitchX(nTemp, n),  m_dMAssmPitchY(nTemp, n), dHeight, assm), *igeomImpl);
01936 
01937 
01938           // position the outer block to match the pins
01939           dX = m_dMAssmPitchX(nTemp, n)/4.0;
01940           dY =  m_dMAssmPitchY(nTemp, n)/4.0;
01941           dZ = (m_dMZAssm(nTemp, 2) + m_dMZAssm(nTemp, 1))/2.0;
01942           //    std::cout << "Move " <<   dMoveX << " " << dMoveY <<std::endl;
01943           IBERRCHK(igeomImpl->moveEnt(assm, dMoveX,dMoveY,dZ), *igeomImpl);
01944 
01945 
01946 
01947           // populate the outer covering array squares
01948           assms[nCount]=assm;
01949         }
01950       }
01951     }
01952 
01953   }
01954 
01955   void AssyGen::CreateOuterCovering ()
01956   // ---------------------------------------------------------------------------
01957   // Function: this function sets the names of the coverings
01958   // Input:    error code
01959   // Output:   none
01960   // ---------------------------------------------------------------------------
01961   {
01962     double xmin, xmax, ymin, ymax, zmin, zmax;
01963     iBase_TagHandle this_tag;
01964     char* tag_name =(char *)"NAME";
01965     std::string sMatName = "";
01966     std::string sMatName0 = "";
01967     std::string sMatName1 = "";
01968 
01969     // get tag handle for 'NAME' tag, already created as iGeom instance is created
01970     IBERRCHK(igeomImpl->getTagHandle(tag_name, this_tag), *igeomImpl);
01971 
01972     iBase_EntityHandle tmp_vol= NULL, tmp_new= NULL;
01973 
01974     // name the innermost outer covering common for both rectangular and hexagonal assembliees
01975     if(m_nDimensions >0){
01976       for (int nTemp1 = 1; nTemp1 <=m_nDuct; nTemp1++){
01977         for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
01978           if(strcmp ( m_szMMAlias(nTemp1, 1).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
01979             sMatName =  m_szAssmMat(p);
01980           }
01981         }
01982 
01983         std::cout << "\ncreated innermost block: " << sMatName << std::endl;
01984 
01985         tmp_vol = assms[(nTemp1 - 1)*m_nDimensions];
01986         IBERRCHK(igeomImpl->setData(tmp_vol, this_tag, sMatName.c_str()), *igeomImpl);
01987 
01988 
01989         Name_Faces( sMatName, tmp_vol, this_tag);
01990       }
01991 
01992       int count =0;//index for edge names
01993       for (int nTemp = 1; nTemp <= m_nDuct; nTemp++){
01994         //  Naming outermost block edges - sidesets in cubit journal file
01995         std::cout << "Naming outermost block edges" << std::endl;
01996         //SimpleArray
01997         std::vector<iBase_EntityHandle> edges;
01998         IBERRCHK(igeomImpl->getEntAdj(assms[nTemp*m_nDimensions -1], iBase_EDGE, edges), *igeomImpl);
01999 
02000 
02001         // get the top corner edges of the outer most covering
02002         std::ostringstream os;
02003         for (int i = 0; i < (int) edges.size(); ++i){
02004           IBERRCHK(igeomImpl->getEntBoundBox(edges[i], xmin, ymin,zmin, xmax, ymax, zmax), *igeomImpl);
02005 
02006 
02007           double dTol = 1e-2; // tolerance for comparing coordinates
02008 
02009           if(abs(zmax - m_dMZAssm(nTemp, 2)) <  dTol){
02010             if(abs(zmax-zmin) < dTol){
02011 
02012               //we have a corner edge - name it
02013               sMatName="side_edge";
02014               ++count;
02015               os << sMatName << count;
02016               sMatName=os.str();
02017               tmp_vol=edges[i];
02018               IBERRCHK(igeomImpl->setData(tmp_vol, this_tag, sMatName.c_str()), *igeomImpl);
02019 
02020               std::cout << "created: " << sMatName << std::endl;
02021               os.str("");
02022               sMatName="";
02023             }
02024           }
02025         }
02026       }
02027       // now subtract the outermost hexes and name them
02028       std::cout << "Subtract outermost hexes and naming them" << std::endl;
02029       int nCount = 0;
02030       for(int nTemp=1; nTemp<=m_nDuct; nTemp++){
02031         for(int n=m_nDimensions; n>1 ; n--){
02032           if(n>1){
02033             ++nCount;
02034             // copy cyl before subtract
02035             IBERRCHK(igeomImpl->copyEnt(assms[(nTemp-1)*m_nDimensions + n-2], tmp_vol), *igeomImpl);
02036 
02037 
02038             // subtract outer most cyl from brick
02039             IBERRCHK(igeomImpl->subtractEnts(assms[(nTemp-1)*m_nDimensions + n-1], tmp_vol, tmp_new), *igeomImpl);
02040 
02041 
02042             assms[(nTemp-1)*m_nDimensions + n-1]=tmp_new;
02043 
02044             // name the vols by searching for the full name of the abbreviated Cell Mat
02045             for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
02046               if(strcmp ( m_szMMAlias(nTemp, n).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02047                 sMatName =  m_szAssmMat(p);
02048               }
02049             }
02050             std::cout << "created: " << sMatName << std::endl;
02051             IBERRCHK(igeomImpl->setData(tmp_new, this_tag, sMatName.c_str()), *igeomImpl);
02052 
02053             Name_Faces( sMatName, tmp_new, this_tag);
02054           }
02055         }
02056       }
02057       std::cout << "\n--------------------------------------------------"<<std::endl;
02058     }
02059 
02060   }
02061 
02062   void AssyGen::Subtract_Pins()
02063   // ---------------------------------------------------------------------------
02064   // Function: subtract the pins from the outer block
02065   // Input:    none
02066   // Output:   none
02067   // ---------------------------------------------------------------------------
02068   {
02069     if (m_nDimensions >0 && in_pins.size()>0){
02070       std::vector<iBase_EntityHandle> copy_inpins(in_pins.size());
02071       std::cout <<"Total number of pins in the model = " <<  in_pins.size()/m_nDuct << std::endl;
02072 
02073       int num_inpins = in_pins.size()/m_nDuct;
02074       CMatrix<iBase_EntityHandle> cp_inpins(m_nDuct, num_inpins);
02075 
02076       for (int k=1; k<=m_nDuct; k++){
02077 
02078         // put all the in pins in a matrix of size duct for subtraction with ducts
02079         for (int i=0; i<num_inpins; i++){
02080           IBERRCHK(igeomImpl->copyEnt(in_pins[ k - 1 + m_nDuct*i], cp_inpins(k,i+1)), *igeomImpl);
02081 
02082         }
02083 
02084         iBase_EntityHandle unite= NULL, tmp_vol = NULL, tmp_new1 = NULL;
02085         tmp_vol = assms[(k-1)*m_nDimensions];
02086 
02087         // subtract the innermost hex from the pins
02088         std::cout << "Duct no.: " << k << " subtracting " << num_inpins << " pins from the duct .. " << std::endl;
02089 
02090         // if there are more than one pins
02091         if(in_pins.size() > 1){
02092           IBERRCHK(igeomImpl->uniteEnts(&cp_inpins(k,1), num_inpins, unite), *igeomImpl);
02093 
02094 
02095           IBERRCHK(igeomImpl->subtractEnts(tmp_vol,unite, tmp_new1), *igeomImpl);
02096 
02097 
02098           tmp_vol = tmp_new1;
02099           unite = NULL;
02100           tmp_new1=NULL;
02101         }
02102         else{ // only one pin in in_pins
02103           IBERRCHK(igeomImpl->subtractEnts(tmp_vol, in_pins[0], tmp_new1), *igeomImpl);
02104 
02105         }
02106       }
02107       std::cout << "\n--------------------------------------------------"<<std::endl;
02108     }
02109     else{
02110       std::cout <<"Nothing to subtract" << std::endl;
02111     }
02112 
02113   }
02114 
02115   void AssyGen::Imprint_Merge()
02116   // ---------------------------------------------------------------------------
02117   // Function: Imprint and Merge
02118   // Input:    none
02119   // Output:   none
02120   // ---------------------------------------------------------------------------
02121   {
02122     // getting all entities for merge and imprint
02123     //SimpleArray
02124     std::vector<iBase_EntityHandle> entities;
02125     IBERRCHK(igeomImpl->getEntities(root_set, iBase_REGION, entities), *igeomImpl);
02126 
02127 
02128     //  now imprint
02129     std::cout << "\n\nImprinting...." << std::endl;
02130     IBERRCHK(igeomImpl->imprintEnts(&entities[0], (int) entities.size()), *igeomImpl);
02131 
02132     std::cout << "\n--------------------------------------------------"<<std::endl;
02133 
02134     //   // merge tolerance
02135     //   double dTol = 1e-4;
02136     //   // now  merge
02137     //   std::cout << "\n\nMerging...." << std::endl;
02138     //   iGeom_mergeEnts(geom, ARRAY_IN(entities), dTol, &err);
02139     //   CHECK("Merge failed.");
02140     //   std::cout <<"merging finished."<< std::endl;
02141     //   std::cout << "\n--------------------------------------------------"<<std::endl;
02142 
02143   }
02144 
02145   void AssyGen::Create2DSurf ()
02146   // ---------------------------------------------------------------------------
02147   // Function: creating planar top surface with zmax
02148   // Input:    error code
02149   // Output:   none
02150   // ---------------------------------------------------------------------------
02151   {
02152     //SimpleArray<iBase_EntityHandle>  all_geom;
02153     //SimpleArray<iBase_EntityHandle> surfs;
02154     int *offset = NULL;
02155     int t=0;
02156     std::cout << "Creating surface; 2D assembly specified..." << std::endl;
02157 
02158     // get all the entities in the model (delete after making a copy of top surface)
02159     std::vector<iBase_EntityHandle> all_geom, surfs;
02160     IBERRCHK(igeomImpl->getEntities(root_set, iBase_REGION, all_geom), *igeomImpl);
02161 
02162 
02163     // get all the surfaces in the model
02164     IBERRCHK(igeomImpl->getArrAdj(&all_geom[0], (int) all_geom.size(), iBase_FACE, surfs, offset), *igeomImpl);
02165 
02166 
02167     double *max_corn = new double [3*surfs.size()];
02168     double *min_corn = new double [3*surfs.size()];
02169     //, min_corn[surfs.size()];
02170     IBERRCHK(igeomImpl->getArrBoundBox(&surfs[0], (int) surfs.size(),iBase_INTERLEAVED, &min_corn[0], &max_corn[0]), *igeomImpl);
02171 
02172 
02173     // find the number of surfaces 't' for array allocation
02174     int nTemp = 1;
02175     double dTol = 1e-3;
02176     double dtop = m_dMZAssm(nTemp, 2);
02177     for (int i = 0; i < (int) surfs.size(); ++i){
02178       if((abs(max_corn[3*i+2] -  dtop) < dTol) && (abs(min_corn[3*i+2] - dtop)<dTol))
02179         t++;
02180     }
02181 
02182     // allocate arrays
02183     SimpleArray<iBase_EntityHandle> max_surfs(t);
02184     SimpleArray<iBase_EntityHandle> new_surfs(t);
02185     t=0;
02186 
02187     // store the max surfaces in max_surfs
02188     for (int i = 0; i < (int) surfs.size(); ++i){
02189 
02190       // locate surfaces for which max and min zcoord is same as maxz coord
02191       if((abs(max_corn[3*i+2] -  dtop) < dTol) && (abs(min_corn[3*i+2] - dtop) < dTol)){
02192         max_surfs[t] = surfs[i];
02193         t++;
02194       }
02195     }
02196 
02197     // make a copy of max_surfs
02198     for(int i = 0; i < max_surfs.size(); ++i){
02199       IBERRCHK(igeomImpl->copyEnt(max_surfs[i], new_surfs[i]), *igeomImpl);
02200 
02201     }
02202 
02203     // delete all the old ents
02204     for(int i=0; i< (int) all_geom.size(); i++){
02205       IBERRCHK(igeomImpl->deleteEnt(all_geom[i]), *igeomImpl);
02206     }
02207     // position the final assembly at the center
02208     // get the assembly on z=0 plane
02209     double zcenter = m_dMZAssm(nTemp, 2)/2.0;//move up
02210     //SimpleArray
02211     std::vector<iBase_EntityHandle> all;
02212 
02213     IBERRCHK(igeomImpl->getEntities(root_set, iBase_REGION, all_geom), *igeomImpl);
02214 
02215 
02216     for(int i=0; i< (int) all.size(); i++){
02217       IBERRCHK(igeomImpl->moveEnt(all[i],0,0,-zcenter), *igeomImpl);
02218 
02219     }
02220     std::cout << "--------------------------------------------------"<<std::endl;
02221 
02222     free(offset);
02223 
02224   }
02225 
02226 
02227   void AssyGen::CreatePinCell( int i, double dX, double dY, double dZ)
02228   //---------------------------------------------------------------------------
02229   //Function: Create pincell i in location dX dY and dZ
02230   //Input:    none
02231   //Output:   none
02232   //---------------------------------------------------------------------------
02233   {
02234     int nRadii=0, nCyl=0, nCells = 0;
02235     double dCylMoveX = 0.0, dCylMoveY = 0.0, dHeightTotal = 0.0;
02236     double dHeight =0.0,dZMove = 0.0, PX = 0.0,PY = 0.0,PZ = 0.0, dP=0.0;
02237     CVector<double> dVCylZPos(2), dVCylXYPos(2), dVStartZ, dVEndZ;;
02238     CVector<std::string> szVMatName, szVMatAlias, szVCellMat;
02239     iBase_EntityHandle cell = NULL, cyl= NULL, tmp_vol= NULL,tmp_vol1= NULL, tmp_new= NULL;
02240 
02241     // name tag handle
02242     iBase_TagHandle this_tag= NULL;
02243     char* tag_name = (char*)"NAME";
02244 
02245     std::string sMatName = "";
02246     std::string sMatName1 = "";
02247 
02248     // get tag handle for 'NAME' tag, already created as iGeom instance is created
02249     IBERRCHK(igeomImpl->getTagHandle(tag_name, this_tag), *igeomImpl);
02250 
02251 
02252 
02253     // get cell material
02254     m_Pincell(i).GetCellMatSize(nCells);
02255     SimpleArray<iBase_EntityHandle> cells(nCells);
02256 
02257     // branch when cells are present
02258     if(nCells > 0){
02259       dVStartZ.SetSize(nCells);
02260       dVEndZ.SetSize(nCells);
02261       szVCellMat.SetSize(nCells);
02262       m_Pincell(i).GetCellMat(dVStartZ, dVEndZ, szVCellMat);
02263 
02264       // get cylinder data
02265       m_Pincell(i).GetNumCyl(nCyl);
02266 
02267       for(int n=1;n<=nCells; n++){
02268 
02269         dHeight = fabs(dVEndZ(n) - dVStartZ(n));
02270 
02271         if(m_szGeomType =="hexagonal"){
02272 
02273           m_Pincell(i).GetPitch(dP, dHeightTotal); // this dHeight is not used in creation
02274           double dSide = dP/(sqrt(3));
02275 
02276           if(nCells >0){
02277             // create prism
02278             IBERRCHK(igeomImpl->createPrism( dHeight, 6, dSide, dSide, cell), *igeomImpl);
02279 
02280           }
02281         }
02282         // if rectangular geometry
02283         if(m_szGeomType =="rectangular"){
02284 
02285           m_Pincell(i).GetPitch(PX, PY, PZ);
02286 
02287           if(nCells >0){
02288             // create brick
02289             IBERRCHK(igeomImpl->createBrick(PX,PY,dHeight,cell), *igeomImpl);
02290 
02291           }
02292         }
02293 
02294         dZMove = (dVEndZ(n)+dVEndZ(n-1))/2.0;
02295 
02296         if(nCells > 0){
02297           // position the brick in assembly
02298           IBERRCHK(igeomImpl->moveEnt( cell, dX, dY, dZMove), *igeomImpl);
02299 
02300           cells[n-1]=cell;
02301 
02302           //search for the full name of the abbreviated Cell Mat and set name
02303           for(int p=1;p<= m_szAssmMatAlias.GetSize();p++){
02304             if(strcmp (szVCellMat(n).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02305               sMatName = m_szAssmMat(p);
02306             }
02307           }
02308 
02309           std::cout << "created: " << sMatName << std::endl;
02310           IBERRCHK(igeomImpl->setData(cell, this_tag, sMatName.c_str()), *igeomImpl);
02311 
02312 
02313           Name_Faces( sMatName, cell, this_tag);
02314         }
02315         // loop and create cylinders
02316         if(nCyl > 0){
02317           m_Pincell(i).GetCylSizes(n, nRadii);
02318           SimpleArray<iBase_EntityHandle> cyls(nRadii);
02319 
02320           //declare variables
02321           CVector<double> dVCylRadii(nRadii);
02322           CVector<std::string> szVMat(nRadii);
02323           CVector<std::string> szVCylMat(nRadii);
02324 
02325           //get values
02326           m_Pincell(i).GetCylRadii(n, dVCylRadii);
02327           m_Pincell(i).GetCylPos(n, dVCylXYPos);
02328           m_Pincell(i).GetCylMat(n, szVCylMat);
02329           m_Pincell(i).GetCylZPos(n, dVCylZPos);
02330           dHeight = dVCylZPos(2)-dVCylZPos(1);
02331 
02332           for (int m=1; m<=nRadii; m++){
02333             IBERRCHK(igeomImpl->createCylinder(dHeight, dVCylRadii(m), dVCylRadii(m), cyl), *igeomImpl);
02334 
02335 
02336             // move their centers and also move to the assembly location  ! Modify if cyl is outside brick
02337             dCylMoveX = dVCylXYPos(1)+dX;
02338             dCylMoveY = dVCylXYPos(2)+dY;
02339             dZMove = (dVCylZPos(1)+dVCylZPos(2))/2.0;
02340             IBERRCHK(igeomImpl->moveEnt(cyl, dCylMoveX,dCylMoveY,dZMove), *igeomImpl);
02341 
02342             ;
02343             cyls[m-1] = cyl;
02344           }
02345 
02346           if(nCells > 0){
02347             // copy cyl before subtract
02348             IBERRCHK(igeomImpl->copyEnt(cyls[nRadii-1], tmp_vol), *igeomImpl);
02349 
02350 
02351             // subtract outer most cyl from brick
02352             IBERRCHK(igeomImpl->subtractEnts( cells[n-1], cyls[nRadii-1], tmp_new), *igeomImpl);
02353 
02354 
02355             // copy the new into the cyl array
02356             cells[n-1] = tmp_new; cell = tmp_new;
02357             cyls[nRadii-1]=tmp_vol;
02358 
02359           }
02360           //set tag on inner most cylinder, search for the full name of the abbreviated Cell Mat
02361           for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
02362             if(strcmp (szVCylMat(1).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02363               sMatName = m_szAssmMat(p);
02364             }
02365           }
02366           tmp_vol1=cyls[0]; //inner most cyl
02367           IBERRCHK(igeomImpl->setData(tmp_vol1, this_tag, sMatName.c_str()), *igeomImpl);
02368 
02369           Name_Faces( sMatName, tmp_vol1, this_tag);
02370 
02371           // other cyl annulus after substraction
02372           for (int b=nRadii; b>1; b--){
02373             IBERRCHK(igeomImpl->copyEnt(cyls[b-2], tmp_vol), *igeomImpl);
02374 
02375 
02376             //subtract tmp vol from the outer most
02377             IBERRCHK(igeomImpl->subtractEnts( cyls[b-1], tmp_vol, tmp_new), *igeomImpl);
02378 
02379 
02380 
02381             // now search for the full name of the abbreviated Cell Mat
02382             for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
02383               if(strcmp (szVCylMat(b).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02384                 sMatName =  m_szAssmMat(p);
02385               }
02386             }
02387             std::cout << "created: " << sMatName << std::endl;
02388             // set the name of the annulus
02389             IBERRCHK(igeomImpl->setData(tmp_new, this_tag, sMatName.c_str()), *igeomImpl);
02390 
02391             Name_Faces( sMatName, tmp_new, this_tag);
02392             // copy the new into the cyl array
02393             cyls[b-1] = tmp_new;
02394           }
02395         }
02396       }
02397     }
02398     // this branch of the routine is responsible for creating cylinders with '0' cells
02399     if(nCells == 0){
02400 
02401       // get cylinder data
02402       m_Pincell(i).GetNumCyl(nCyl);
02403       nCells = nCyl;
02404 
02405 
02406       for(int n=1;n<=nCells; n++){
02407 
02408 
02409         if(m_szGeomType =="hexagonal"){
02410 
02411           m_Pincell(i).GetPitch(dP, dHeightTotal); // this dHeight is not used in creation
02412         }
02413         // if rectangular geometry
02414         if(m_szGeomType =="rectangular"){
02415 
02416           m_Pincell(i).GetPitch(PX, PY, PZ);
02417         }
02418 
02419         // loop and create cylinders
02420         if(nCyl > 0){
02421           m_Pincell(i).GetCylSizes(n, nRadii);
02422           SimpleArray<iBase_EntityHandle> cyls(nRadii);
02423 
02424           //declare variables
02425           CVector<double> dVCylRadii(nRadii);
02426           CVector<std::string> szVMat(nRadii);
02427           CVector<std::string> szVCylMat(nRadii);
02428 
02429           //get values
02430           m_Pincell(i).GetCylRadii(n, dVCylRadii);
02431           m_Pincell(i).GetCylPos(n, dVCylXYPos);
02432           m_Pincell(i).GetCylMat(n, szVCylMat);
02433           m_Pincell(i).GetCylZPos(n, dVCylZPos);
02434 
02435           dHeight = dVCylZPos(2)-dVCylZPos(1);
02436 
02437           for (int m=1; m<=nRadii; m++){
02438             IBERRCHK(igeomImpl->createCylinder(dHeight, dVCylRadii(m), dVCylRadii(m), cyl), *igeomImpl);
02439 
02440 
02441             // move their centers and also move to the assembly location  ! Modify if cyl is outside brick
02442             dCylMoveX = dVCylXYPos(1)+dX;
02443             dCylMoveY = dVCylXYPos(2)+dY;
02444             dZMove = (dVCylZPos(1)+dVCylZPos(2))/2.0;
02445             IBERRCHK(igeomImpl->moveEnt( cyl, dCylMoveX,dCylMoveY,dZMove), *igeomImpl);
02446 
02447 
02448             cyls[m-1] = cyl;
02449           }
02450 
02451           //set tag on inner most cylinder, search for the full name of the abbreviated Cell Mat
02452           for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
02453             if(strcmp (szVCylMat(1).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02454               sMatName = m_szAssmMat(p);
02455             }
02456           }
02457           std::cout << "created: " << sMatName << std::endl;
02458           tmp_vol1=cyls[0]; //inner most cyl
02459 
02460           IBERRCHK(igeomImpl->setData(tmp_vol1, this_tag, sMatName.c_str()), *igeomImpl);
02461 
02462 
02463           Name_Faces( sMatName, tmp_vol1, this_tag);
02464 
02465           // other cyl annulus after substraction
02466           for (int b=nRadii; b>1; b--){
02467             IBERRCHK(igeomImpl->copyEnt(cyls[b-2], tmp_vol), *igeomImpl);
02468 
02469 
02470             //subtract tmp vol from the outer most
02471             IBERRCHK(igeomImpl->subtractEnts( cyls[b-1], tmp_vol, tmp_new), *igeomImpl);
02472 
02473 
02474 
02475             // now search for the full name of the abbreviated Cell Mat
02476             for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
02477               if(strcmp (szVCylMat(b).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02478                 sMatName =  m_szAssmMat(p);
02479               }
02480             }
02481             std::cout << "created: " << sMatName << std::endl;
02482             // set the name of the annulus
02483             IBERRCHK(igeomImpl->setData(tmp_new, this_tag, sMatName.c_str()), *igeomImpl);
02484 
02485             Name_Faces( sMatName, tmp_new, this_tag);
02486 
02487             // copy the new into the cyl array
02488             cyls[b-1] = tmp_new;
02489           }
02490         }
02491       }
02492     }
02493 
02494   }
02495 
02496 
02497   void AssyGen::CreatePinCell_Intersect( int i, double dX, double dY, double dZ)
02498   //---------------------------------------------------------------------------
02499   //Function: Create pincell i in location dX dY and dZ
02500   //Input:    none
02501   //Output:   none
02502   //---------------------------------------------------------------------------
02503   {
02504     int nRadii=0, nCyl=0, nCells = 0;
02505     double dCylMoveX = 0.0, dCylMoveY = 0.0, dHeightTotal = 0.0;
02506     double dHeight =0.0,dZMove = 0.0, PX = 0.0,PY = 0.0,PZ = 0.0, dP=0.0;
02507     CVector<double> dVCylZPos(2), dVCylXYPos(2), dVEndZ, dVStartZ;
02508     CVector<std::string> szVMatName, szVMatAlias, szVCellMat;
02509     iBase_EntityHandle cell = NULL, cyl= NULL, tmp_vol1= NULL, tmp_new= NULL, cell_copy = NULL, intersec = NULL;
02510 
02511     // name tag handle
02512     iBase_TagHandle this_tag= NULL;
02513     char* tag_name = (char*)"NAME";
02514 
02515     std::string sMatName = "";
02516     std::string sMatName0 = "";
02517     std::string sMatName1 = "";
02518 
02519     // get tag handle for 'NAME' tag, already created as iGeom instance is created
02520     IBERRCHK(igeomImpl->getTagHandle(tag_name, this_tag), *igeomImpl);
02521 
02522 
02523     // get cell material
02524     m_Pincell(i).GetCellMatSize(nCells);
02525     m_Pincell(i).GetNumCyl(nCyl);
02526     SimpleArray<iBase_EntityHandle> cells(nCells);
02527 
02528     // branch when cells are present
02529     if(nCells > 0){
02530       dVStartZ.SetSize(nCells);
02531       dVEndZ.SetSize(nCells);
02532       szVCellMat.SetSize(nCells);
02533       m_Pincell(i).GetCellMat(dVStartZ, dVEndZ, szVCellMat);
02534 
02535       // get cylinder data
02536       m_Pincell(i).GetNumCyl(nCyl);
02537 
02538       for(int n=1;n<=nCells; n++){
02539 
02540         dHeight = dVEndZ(n) - dVStartZ(n);
02541 
02542         if(m_szGeomType =="hexagonal"){
02543 
02544           m_Pincell(i).GetPitch(dP, dHeightTotal); // this dHeight is not used in creation
02545 
02546           double dSide = dP/(sqrt(3));
02547 
02548           if(nCells >0){
02549             // create prism
02550             IBERRCHK(igeomImpl->createPrism( dHeight, 6, dSide, dSide, cell), *igeomImpl);
02551 
02552           }
02553         }
02554         // if rectangular geometry
02555         if(m_szGeomType =="rectangular"){
02556 
02557           m_Pincell(i).GetPitch(PX, PY, PZ);
02558 
02559           if(nCells >0){
02560             // create brick
02561             IBERRCHK(igeomImpl->createBrick(PX,PY,dHeight, cell), *igeomImpl);
02562 
02563           }
02564         }
02565 
02566         dZMove = (dVEndZ(n)+dVEndZ(n-1))/2.0;
02567 
02568         if(nCells > 0){
02569           // position the brick in assembly
02570           IBERRCHK(igeomImpl->moveEnt(  cell, dX, dY, dZMove), *igeomImpl);
02571 
02572 
02573           cells[n-1]=cell;
02574         }
02575         // loop and create cylinders
02576         if(nCyl > 0){
02577           m_Pincell(i).GetCylSizes(n, nRadii);
02578           SimpleArray<iBase_EntityHandle> cyls(nRadii);
02579           SimpleArray<iBase_EntityHandle> cell_copys(nRadii);
02580           SimpleArray<iBase_EntityHandle> intersec_main(nRadii);
02581           iBase_EntityHandle  tmp_intersec;
02582           //declare variables
02583           CVector<double> dVCylRadii(nRadii);
02584           CVector<std::string> szVMat(nRadii);
02585           CVector<std::string> szVCylMat(nRadii);
02586 
02587           //get values
02588           m_Pincell(i).GetCylRadii(n, dVCylRadii);
02589           m_Pincell(i).GetCylPos(n, dVCylXYPos);
02590           m_Pincell(i).GetCylMat(n, szVCylMat);
02591           m_Pincell(i).GetCylZPos(n, dVCylZPos);
02592           dHeight = dVCylZPos(2)-dVCylZPos(1);
02593 
02594           for (int m=1; m<=nRadii; m++){
02595             IBERRCHK(igeomImpl->createCylinder(dHeight, dVCylRadii(m), dVCylRadii(m), cyl), *igeomImpl);
02596 
02597 
02598             // move their centers and also move to the assembly location  ! Modify if cyl is outside brick
02599             dCylMoveX = dVCylXYPos(1)+dX;
02600             dCylMoveY = dVCylXYPos(2)+dY;
02601             dZMove = (dVCylZPos(1)+dVCylZPos(2))/2.0;
02602             IBERRCHK(igeomImpl->moveEnt( cyl, dCylMoveX,dCylMoveY,dZMove), *igeomImpl);
02603 
02604 
02605             cyls[m-1] = cyl;
02606 
02607 
02608             //copy cell nRadii  times for intersection with cylinders
02609             IBERRCHK(igeomImpl->copyEnt(cells[n-1], cell_copy), *igeomImpl);
02610 
02611 
02612             cell_copys[m-1] = cell_copy;
02613             IBERRCHK(igeomImpl->intersectEnts(cell_copys[m-1], cyls[m-1], intersec), *igeomImpl);
02614 
02615 
02616             intersec_main[m-1] = intersec;
02617             intersec = NULL;
02618           }
02619 
02620           //set tag on inner most cylinder, search for the full name of the abbreviated Cell Mat
02621           tmp_vol1=intersec_main[0];
02622           for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
02623             if(strcmp (szVCylMat(1).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02624               sMatName = m_szAssmMat(p);
02625             }
02626           }
02627           IBERRCHK(igeomImpl->setData(tmp_vol1, this_tag, sMatName.c_str()), *igeomImpl);
02628 
02629 
02630           Name_Faces( sMatName, tmp_vol1, this_tag);
02631 
02632           // copy the outermost cyl
02633           IBERRCHK(igeomImpl->copyEnt(intersec_main[nRadii-1], tmp_intersec), *igeomImpl);
02634 
02635 
02636           // subtract the outermost cyl from the cell
02637           IBERRCHK(igeomImpl->subtractEnts( cells[n-1], tmp_intersec, tmp_new), *igeomImpl);
02638 
02639 
02640           // now search for the full name of the abbreviated Cell Mat
02641           for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
02642             if(strcmp (szVCellMat(n).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02643               sMatName =  m_szAssmMat(p);
02644             }
02645           }
02646           std::cout << "created: " << sMatName << std::endl;
02647           // set the name of the annulus
02648           IBERRCHK(igeomImpl->setData(tmp_new, this_tag, sMatName.c_str()), *igeomImpl);
02649 
02650           Name_Faces( sMatName, tmp_new, this_tag);
02651 
02652           for (int b=nRadii; b>1; b--){
02653             IBERRCHK(igeomImpl->copyEnt(intersec_main[b-2], tmp_intersec), *igeomImpl);
02654 
02655 
02656             //subtract tmp vol from the outer most
02657             IBERRCHK(igeomImpl->subtractEnts( intersec_main[b-1], tmp_intersec, tmp_new), *igeomImpl);
02658 
02659 
02660 
02661             // now search for the full name of the abbreviated Cell Mat
02662             for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
02663               if(strcmp (szVCylMat(b).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02664                 sMatName =  m_szAssmMat(p);
02665               }
02666             }
02667             std::cout << "created: " << sMatName << std::endl;
02668             // set the name of the annulus
02669             IBERRCHK(igeomImpl->setData(tmp_new, this_tag, sMatName.c_str()), *igeomImpl);
02670 
02671             Name_Faces( sMatName, tmp_new, this_tag);
02672 
02673             // copy the new into the cyl array
02674             cyls[b-1] = tmp_new;
02675 
02676           }
02677         }
02678       }
02679     }
02680     // this branch of the routine is responsible for creating cylinders with '0' cells
02681     if(nCells == 0){
02682 
02683       // get cylinder data
02684       m_Pincell(i).GetNumCyl(nCyl);
02685       nCells = nCyl;
02686       cells.resize(nCells);
02687 
02688       for(int n=1;n<=nCells; n++){
02689 
02690         // get some cylinder parameters to create the cell material for intersection
02691         m_Pincell(i).GetCylZPos(n, dVCylZPos);
02692         dHeight = dVCylZPos(2)-dVCylZPos(1);
02693         dZMove = (dVCylZPos(1)+dVCylZPos(2))/2.0;
02694 
02695         if(m_szGeomType =="hexagonal"){
02696 
02697           m_Pincell(i).GetPitch(dP, dHeightTotal); // this dHeight is not used in creation
02698           double dSide = dP/(sqrt(3));
02699           IBERRCHK(igeomImpl->createPrism(dHeight, 6, dSide, dSide, cell), *igeomImpl);
02700 
02701 
02702         }
02703         // if rectangular geometry
02704         if(m_szGeomType =="rectangular"){
02705 
02706           m_Pincell(i).GetPitch(PX, PY, PZ);
02707           // create brick
02708           IBERRCHK(igeomImpl->createBrick(PX,PY,dHeight,cell), *igeomImpl);
02709 
02710 
02711         }
02712         IBERRCHK(igeomImpl->moveEnt( cell, dX, dY, dZMove), *igeomImpl);
02713 
02714 
02715 
02716         cells[n-1]=cell;
02717         // loop and create cylinders
02718         if(nCyl > 0){
02719           m_Pincell(i).GetCylSizes(n, nRadii);
02720 
02721           //declare variables
02722           SimpleArray<iBase_EntityHandle> cyls(nRadii), cell_copys(nRadii), intersec_main(nRadii), intersec_copy(nRadii);
02723           iBase_EntityHandle  tmp_intersec;
02724           CVector<double> dVCylRadii(nRadii);
02725           CVector<std::string> szVMat(nRadii), szVCylMat(nRadii);
02726 
02727           //get values
02728           m_Pincell(i).GetCylRadii(n, dVCylRadii);
02729           m_Pincell(i).GetCylPos(n, dVCylXYPos);
02730           m_Pincell(i).GetCylMat(n, szVCylMat);
02731 
02732           for (int m=1; m<=nRadii; m++){
02733             IBERRCHK(igeomImpl->createCylinder(dHeight, dVCylRadii(m), dVCylRadii(m), cyl), *igeomImpl);
02734 
02735 
02736             // move their centers and also move to the assembly location  ! Modify if cyl is outside brick
02737             dCylMoveX = dVCylXYPos(1)+dX;
02738             dCylMoveY = dVCylXYPos(2)+dY;
02739 
02740             IBERRCHK(igeomImpl->moveEnt( cyl, dCylMoveX,dCylMoveY,dZMove), *igeomImpl);
02741 
02742 
02743             cyls[m-1] = cyl;
02744 
02745             //copy cell nRadii  times for intersection with cylinders
02746             IBERRCHK(igeomImpl->copyEnt(cells[n-1], cell_copy), *igeomImpl);
02747 
02748             //    cell_copys[m-1] = cell_copy;
02749             IBERRCHK(igeomImpl->intersectEnts(cell_copy, cyls[m-1], intersec), *igeomImpl);
02750 
02751             intersec_main[m-1] = intersec;
02752             intersec = NULL;
02753           }
02754 
02755           //set tag on inner most cylinder, search for the full name of the abbreviated Cell Mat
02756           tmp_vol1=intersec_main[0];
02757           for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
02758             if(strcmp (szVCylMat(1).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02759               sMatName = m_szAssmMat(p);
02760             }
02761           }
02762           IBERRCHK(igeomImpl->setData(tmp_vol1, this_tag, sMatName.c_str()), *igeomImpl);
02763 
02764 
02765           Name_Faces( sMatName, tmp_vol1, this_tag);
02766 
02767           // delete the cell as this is the case when no. cell material is specified
02768           IBERRCHK(igeomImpl->deleteEnt(cells[n-1]), *igeomImpl);
02769 
02770 
02771 
02772           // other cyl annulus after substraction
02773           for (int b=nRadii; b>1; b--){
02774             IBERRCHK(igeomImpl->copyEnt(intersec_main[b-2], tmp_intersec), *igeomImpl);
02775 
02776             //subtract tmp vol from the outer most
02777             IBERRCHK(igeomImpl->subtractEnts( intersec_main[b-1], tmp_intersec, tmp_new), *igeomImpl);
02778 
02779 
02780 
02781             // now search for the full name of the abbreviated Cell Mat
02782             for(int p=1;p<=m_szAssmMatAlias.GetSize();p++){
02783               if(strcmp (szVCylMat(b).c_str(), m_szAssmMatAlias(p).c_str()) == 0){
02784                 sMatName =  m_szAssmMat(p);
02785               }
02786             }
02787             std::cout << "created: " << sMatName << std::endl;
02788             // set the name of the annulus
02789             IBERRCHK(igeomImpl->setData(tmp_new, this_tag, sMatName.c_str()), *igeomImpl);
02790 
02791             Name_Faces( sMatName, tmp_new, this_tag);
02792 
02793             // copy the new into the cyl array
02794             cyls[b-1] = tmp_new;
02795 
02796           }
02797         }
02798 
02799       }
02800     }
02801 
02802   }
02803 
02804 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines