MeshKit  1.0
EBMesher.cpp
Go to the documentation of this file.
00001 #include "meshkit/EBMesher.hpp"
00002 #include "meshkit/SCDMesh.hpp"
00003 #include "meshkit/MKCore.hpp"
00004 #include "meshkit/ModelEnt.hpp"
00005 #include "meshkit/SizingFunction.hpp"
00006 #include "moab/ReadUtilIface.hpp"
00007 #include "meshkit/iGeom.hpp"
00008 #include "meshkit/RegisterMeshOp.hpp"
00009 
00010 #include "moab/Core.hpp"
00011 #include "moab/Range.hpp"
00012 #include "moab/CartVect.hpp"
00013 #include "moab/EntityType.hpp"
00014 #include "moab/HomXform.hpp"
00015 #include "moab/GeomTopoTool.hpp"
00016 
00017 #include <time.h>
00018 
00019 #ifdef HAVE_CGM
00020 #include "CubitDefines.h"
00021 #include "GeometryQueryTool.hpp"
00022 #include "RefFace.hpp"
00023 #endif
00024 
00025 #define PI 3.14159265
00026 #define ERRORRF(a) {if (iBase_SUCCESS != err) {std::cerr << a << std::endl; return false;}}
00027 
00028 double rayDir[3][3] = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}};
00029 
00030 const bool debug_ebmesh = false;
00031 inline bool less_intersect(const IntersectDist d1, const IntersectDist d2) {
00032   return d1.distance < d2.distance;
00033 }
00034 
00035 inline bool equal_intersect(const IntersectDist d1, const IntersectDist d2) {
00036   return std::abs(d1.distance - d2.distance) < 10e-7;
00037 }
00038 
00039 namespace MeshKit 
00040 {
00041 // static registration of this  mesh scheme
00042 moab::EntityType EBMesher_tps[] = {moab::MBVERTEX, moab::MBHEX, moab::MBMAXTYPE};
00043 const moab::EntityType* EBMesher::output_types()
00044   { return EBMesher_tps; }
00045 
00046 EBMesher::EBMesher(MKCore *mkcore, const MEntVector &me_vec,
00047                    double size, bool use_geom, int add_layer)
00048   : MeshScheme(mkcore, me_vec), m_nAddLayer(add_layer),  m_dInputSize(size), m_bUseGeom(use_geom)
00049 {
00050   m_mesh = mkcore->imesh_instance()->instance();
00051   m_hRootSet = mkcore->imesh_instance()->getRootSet();
00052   m_nStatus = OUTSIDE;
00053   m_iStartHex = 0;
00054   m_hObbTree = NULL;
00055   m_hTreeRoot = 0;
00056 
00057   // create tags with MOAB for dense tags
00058   int outside = 1;
00059   const void *out = &outside;
00060   m_elemStatusTag = get_tag("ELEM_STATUS_TAG", 1,
00061                             moab::MB_TAG_DENSE, moab::MB_TYPE_INTEGER, out);
00062   
00063   int length = 1;
00064   const void *leng = &length;
00065   m_edgeCutFracLengthTag = get_tag("EDGE_CUT_FRAC_LENGTH_TAG", // # of fractions
00066                                    3,
00067                                    //moab::MB_TAG_SPARSE, // using dense for hdf5 exporting performance
00068                                    moab::MB_TAG_DENSE,
00069                                    moab::MB_TYPE_INTEGER, leng);
00070 
00071   m_edgeCutFracTag = get_various_length_tag("EDGE_CUT_FRAC_TAG",
00072                                             //moab::MB_TAG_SPARSE,
00073                                             moab::MB_TAG_DENSE,
00074                                             moab::MB_TYPE_DOUBLE);
00075 
00076   int m_id = 1;
00077   const void *id = &m_id;
00078   m_volFracLengthTag = get_tag("VOL_FRAC_LENGTH_TAG", 1,
00079                                moab::MB_TAG_SPARSE, moab::MB_TYPE_INTEGER, id);
00080   
00081   m_volFracHandleTag = get_various_length_tag("VOL_FRAC_HANDLE_TAG",
00082                                           moab::MB_TAG_SPARSE, moab::MB_TYPE_INTEGER);
00083 
00084   m_volFracTag = get_various_length_tag("VOL_FRAC_TAG",
00085                                         moab::MB_TAG_SPARSE, moab::MB_TYPE_DOUBLE);
00086   
00087   // set bounding box size tag
00088   double bb_default[6] = { 0., 0., 0., 0., 0., 0. };
00089   m_bbTag = get_tag("BOUNDING_BOX_SIZE", 6,
00090                     moab::MB_TAG_SPARSE, moab::MB_TYPE_DOUBLE, bb_default);
00091   
00092   m_GeomTopoTool = new moab::GeomTopoTool( moab_instance() );
00093 }
00094 
00095 EBMesher::~EBMesher()
00096 {
00097   delete m_GeomTopoTool;
00098 }
00099 
00100 bool EBMesher::can_mesh(ModelEnt *me) 
00101 {
00102   if (me->dimension() == 3) return true;
00103   else return false;
00104 }
00105     
00106 bool EBMesher::add_modelent(ModelEnt *model_ent) 
00107 {
00108     // make sure this represents geometric volumes
00109   if (3 != model_ent->dimension()) 
00110     throw Error(MK_WRONG_DIMENSION, "Wrong dimension entity added to EBMesher.");
00111 
00112   return MeshOp::add_modelent(model_ent);
00113 }
00114 
00115 MeshOp *EBMesher::get_scd_mesher()
00116 {
00117   MeshOpProxy* proxy = NULL;
00118   proxy = MKCore::meshop_proxy("SCDMesh");
00119   if (proxy == NULL) throw Error(MK_FAILURE, "Couldn't find a MeshOp capable of producing the given mesh type.");
00120   return mk_core()->construct_meshop(proxy);
00121 }
00122 
00123 void EBMesher::setup_this()
00124 {
00125   MeshOp *scd_mesher = NULL;
00126   std::vector<MeshOp*> meshops;
00127   bool inserted = false;
00128 
00129   for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
00130     ModelEnt *me = (*sit).first;
00131     bool found = false;
00132     if (!scd_mesher) scd_mesher = get_scd_mesher();
00133     meshops.clear();
00134     me->get_meshops(meshops);
00135     int n_meshops = meshops.size();
00136 
00137     if (n_meshops > 0) {
00138       for (int i = 0; i < n_meshops; i++) {
00139         if (meshops[i] == scd_mesher) {
00140           found = true;
00141           break;
00142         }
00143       }
00144     }
00145 
00146     if (!found) { // if no scd meshop
00147       // add this entity to it, and if first, make sure it's added upstream
00148       scd_mesher->add_modelent(me);
00149       if (!inserted) {
00150         mk_core()->insert_node(scd_mesher, this);
00151         inserted = true;
00152       }
00153     }
00154 
00155     // no need to traverse all geometry for using whole geometry
00156     if (m_bUseWholeGeom) break;
00157   }
00158 
00159   SCDMesh *sm = (SCDMesh*) scd_mesher;
00160   sm->set_interface_scheme(SCDMesh::full);
00161   sm->set_grid_scheme(SCDMesh::cfMesh);
00162   
00163   sm->set_axis_scheme(SCDMesh::cartesian);
00164   sm->set_box_increase_ratio(m_boxIncrease); // add some extra layer to box
00165 
00166   if (m_bUseWholeGeom) sm->set_geometry_scheme(SCDMesh::all);
00167   else sm->set_geometry_scheme(SCDMesh::individual);
00168 
00169   // use mesh based geometry in SCDMesh and make tree to get box dimension
00170   if (m_bUseMeshGeom) {
00171     sm->use_mesh_geometry(true);
00172     sm->set_cart_box_min_max(m_minCoord, m_maxCoord, m_boxIncrease);
00173   }
00174 
00175   // set # of intervals for 3 directions
00176   std::vector<int> fine_i (m_nDiv[0], 1);
00177   sm->set_coarse_i_grid(m_nDiv[0]);
00178   sm->set_fine_i_grid(fine_i);
00179   std::vector<int> fine_j (m_nDiv[1], 1);
00180   sm->set_coarse_j_grid(m_nDiv[1]);
00181   sm->set_fine_j_grid(fine_j);
00182   std::vector<int> fine_k (m_nDiv[2], 1);
00183   sm->set_coarse_k_grid(m_nDiv[2]);
00184   sm->set_fine_k_grid(fine_k);
00185 }
00186 
00187 void EBMesher::execute_this() 
00188 {
00189   iBase_ErrorType  err;
00190   GraphNode *scd_node = other_node(in_arcs());
00191 
00192   double obb_time = .0;
00193   double intersection_time = .0;
00194   double set_info_time = .0;
00195   time_t time1, time2, time3, time4;
00196   unsigned long long mem1, mem3, mem4;
00197   moab_instance()->estimated_memory_use(0, 0, 0, &mem1);
00198   moab::ErrorCode rval;
00199 
00200   if (debug_ebmesh) {
00201     rval = mk_core()->moab_instance()->write_mesh("input.vtk");
00202     MBERRCHK(rval, mk_core()->moab_instance());
00203   }
00204 
00205   time(&time1);
00206   if (m_bUseGeom && !m_bUseMeshGeom) { // construct obb tree for geometry surfaces and volumes by GeomTopoTool
00207     rval = m_GeomTopoTool->construct_obb_trees(m_bUseWholeGeom);
00208     MBERRCHK(rval, mk_core()->moab_instance());
00209   }
00210   time(&time2);
00211   obb_time += difftime(time2, time1);
00212 
00213   for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
00214     // 1. set division
00215     double box_min_max[6];
00216     if (m_bUseWholeGeom) {
00217       static_cast<SCDMesh*> (scd_node)->get_box_dimension(box_min_max, &box_min_max[3]);
00218     }
00219     else {
00220       err = mk_core()->imesh_instance()->getData(reinterpret_cast< iBase_EntityHandle >
00221                                                  (sit ->first->mesh_handle()), m_bbTag, &box_min_max[0]);
00222       IBERRCHK(err, "Failed to get hexes.\n");
00223     }
00224     set_division(box_min_max, &box_min_max[3]);
00225 
00226     // 2. set or construct obb tree
00227     time(&time1);
00228     set_tree_root(sit->first);
00229     time(&time2);
00230     obb_time += difftime(time2, time1);
00231 
00232     // 3. find intersected geometry surfaces by rays
00233     find_intersections();
00234     time(&time3);
00235     intersection_time += difftime(time3, time2);
00236     moab_instance()->estimated_memory_use(0, 0, 0, &mem3);
00237     
00238     // 4. set hex status and boundary hex cut fraction info
00239     set_tag_info();
00240     time(&time4);
00241     set_info_time += difftime(time4, time3);
00242     moab_instance()->estimated_memory_use(0, 0, 0, &mem4);
00243 
00244     if (m_bUseWholeGeom) break;
00245   }
00246   if (debug_ebmesh) {
00247     std::cout << "OBB_tree_construct_time: " << obb_time
00248               << ", intersection_time: " << intersection_time
00249               << ", set_info_time: " << set_info_time << std::endl;
00250     std::cout << "start_memory: " << mem1
00251               //<< ", OBB_tree_construct_moemory: " << mem2
00252               << ", intersection_memory: " << mem3
00253               << ", set_info_memory: " << mem4
00254               << std::endl;
00255   }
00256 }
00257 
00258 void EBMesher::set_num_interval(int* n_interval)
00259 {
00260   for (int i = 0; i < 3; i++) {
00261     m_nDiv[i] = n_interval[i];
00262   }
00263 }
00264 
00265 void EBMesher::set_tree_root(ModelEnt* me)
00266 {
00267   moab::ErrorCode rval;
00268   if (m_bUseGeom) {
00269     if (m_hObbTree == NULL) {
00270       m_hObbTree = m_GeomTopoTool->obb_tree();
00271     }
00272     if (m_bUseWholeGeom) {
00273       if (!m_bUseMeshGeom) m_hTreeRoot = m_GeomTopoTool->get_one_vol_root();
00274     }
00275     else {
00276       rval = m_GeomTopoTool->get_root(me->mesh_handle(), m_hTreeRoot);
00277       MBERRCHK(rval, mk_core()->moab_instance());
00278     }
00279   }
00280   else { // facet data input case
00281     moab::EntityHandle meshset;
00282     if (m_bUseWholeGeom) meshset = 0;
00283     else meshset = me->mesh_handle();
00284 
00285     // get triangles
00286     moab::Range tris;
00287     rval = moab_instance()->get_entities_by_dimension(meshset, 2, tris);
00288     MBERRCHK(rval, mk_core()->moab_instance());
00289     
00290     // make tree
00291     if (m_hObbTree == NULL) {
00292       m_hObbTree = new moab::OrientedBoxTreeTool( moab_instance() );
00293     }
00294     rval = m_hObbTree->build(tris, m_hTreeRoot);
00295     MBERRCHK(rval, mk_core()->moab_instance());
00296   }
00297 }
00298 
00299 void EBMesher::find_intersections()
00300 {
00301   std::cout << "starting to find intersections." << std::endl;
00302   int i;
00303   //m_vnHexStatus.resize(m_nHex, 1); // initialize all hex as outside
00304   m_vnHexStatus.resize(m_nHex, 0); // initialize all hex as inside
00305 
00306   // fire rays to 3 directions
00307   for (i = 0; i < 3; i++) {
00308     //m_vnEdgeStatus[i].resize(m_nDiv[i]*m_nDiv[(i + 1)%3]*m_nDiv[(i + 2)%3], OUTSIDE);
00309     m_vnEdgeStatus[i].resize(m_nDiv[i]*m_nDiv[(i + 1)%3]*m_nDiv[(i + 2)%3], INSIDE);
00310     fire_rays(i);
00311   }
00312   std::cout << "finished to find intersections." << std::endl;
00313 }
00314 
00315 void EBMesher::export_mesh(const char* file_name, bool separate)
00316 { 
00317   // get all hexes
00318   time_t time1, time2, time3;
00319   time(&time1);
00320   int i;
00321   iBase_ErrorType err;
00322   std::vector<iBase_EntityHandle> hex_handles;
00323   err = mk_core()->imesh_instance()->getEntities(m_hRootSet, iBase_REGION,
00324                                                  iMesh_HEXAHEDRON, hex_handles);
00325   IBERRCHK(err, "Failed to get hexes.\n");
00326   /*
00327   int hex_size = hex_handles.size();
00328   int* hex_status = new int[hex_size];
00329   int* hex_status = NULL;
00330   std::vector<int> hex_status(hex_size);
00331   err = mk_core()->imesh_instance()->getIntArrData(&hex_handles[0], hex_size,
00332   m_elemStatusTag, &hex_status[0]);*/
00333 
00334   int error;
00335   int hex_size = hex_handles.size();
00336   int* hex_status = new int[hex_size];
00337   int hex_status_alloc = sizeof(int)*hex_size;
00338   int hex_status_size = 0;
00339   iMesh_getIntArrData(m_mesh, &hex_handles[0],
00340                       hex_size, m_elemStatusTag,
00341                       &hex_status, &hex_status_alloc,
00342                       &hex_status_size, &error);
00343   IBERRCHK(error, "Failed to get hex status.\n");
00344   time(&time2);
00345 
00346   if (separate) {
00347     int n_inside_hex = 0;
00348     int n_outside_hex = 0;
00349     int n_boundary_hex = 0;
00350     int hex_stat;
00351     (void) hex_stat;
00352     std::vector<iBase_EntityHandle> insideHex, outsideHex, bndrHex;
00353     for (i = 0; i < hex_size; i++) {
00354       if (hex_status[i] == 0) {
00355         n_inside_hex++;
00356         hex_stat = 0;
00357         insideHex.push_back(hex_handles[i]);
00358       }
00359       else if (hex_status[i] == 1) {
00360         n_outside_hex++;
00361         hex_stat = 1;
00362         outsideHex.push_back(hex_handles[i]);
00363       }
00364       else if (hex_status[i] == 2) {
00365         n_boundary_hex++;
00366         hex_stat = 2;
00367         bndrHex.push_back(hex_handles[i]);
00368       }
00369       else {
00370         throw Error(MK_FAILURE, "Hex element status should be inside/outside/boundary.");
00371       }
00372     }
00373     
00374     std::cout << "# of exported inside hex:" << n_inside_hex
00375               << ", # of exported outside hex:" << n_outside_hex
00376               << ", # of exported boundary hex:" << n_boundary_hex
00377               << ", geom vol:"
00378               << n_inside_hex*m_dIntervalSize[0]*m_dIntervalSize[1]*m_dIntervalSize[2]
00379               << ", vox vol:" << hex_size*m_dIntervalSize[0]*m_dIntervalSize[1]*m_dIntervalSize[2]
00380               << std::endl;
00381     time(&time3);
00382 
00383     // save inside/outside/boundary elements separately
00384     if (n_inside_hex > 0) {
00385       write_mesh(file_name, 0, &insideHex[0], n_inside_hex);
00386     }
00387     if (n_boundary_hex > 0) {
00388       write_mesh(file_name, 2, &bndrHex[0], n_boundary_hex);
00389     }
00390     if (n_outside_hex > 0) {
00391       write_mesh(file_name, 1, &outsideHex[0], n_outside_hex);
00392     }
00393 
00394     if (debug_ebmesh) {
00395       std::cout << "hex_handle_get_time: "
00396                 << difftime(time2, time1)
00397                 << ", separate_write_time: "
00398                 << difftime(time3, time2)
00399                 << std::endl;
00400     }
00401   }
00402   else {
00403     write_mesh(file_name, 3, &hex_handles[0], hex_size);
00404     time3 = clock();
00405     if (debug_ebmesh) {
00406       std::cout << "hex_handle_get_time: "
00407                 << difftime(time2, time1)
00408                 << ", actual_write_time: "
00409                 << difftime(time3, time2)
00410                 << std::endl;
00411     }
00412   }
00413   delete [] hex_status;
00414 }
00415 
00416 void EBMesher::write_mesh(const char* file_name, int type,
00417                           iBase_EntityHandle* handles, int& n_elem)
00418 {
00419   time_t time1, time2, time3;
00420   time(&time1);
00421   int is_list = 1;
00422   moab::ErrorCode rval;
00423   iBase_EntitySetHandle set;
00424   iBase_ErrorType err;
00425 
00426   err = mk_core()->imesh_instance()->createEntSet(is_list, set);
00427   IBERRCHK(err, "Couldn't create set.");
00428 
00429   err = mk_core()->imesh_instance()->addEntArrToSet(handles, n_elem, set);
00430   IBERRCHK(err, "Couldn't add hexes to set.");
00431   time(&time2);
00432 
00433   std::string out_name;
00434   std::stringstream ss;
00435   if (type == 0) ss << "inside_";
00436   else if (type == 1) ss << "outside_";
00437   else if (type == 2) ss << "boundary_";
00438   ss << file_name;
00439   ss >> out_name;
00440   
00441   rval = moab_instance()->write_mesh(out_name.c_str(),
00442                                      (const moab::EntityHandle*) &set, 1);
00443   MBERRCHK(rval, mk_core()->moab_instance());
00444   
00445   std::cout << "Elements are exported." << std::endl;
00446   time(&time3);
00447 
00448   if (debug_ebmesh) {
00449     std::cout << "set_creation_time: "
00450               << difftime(time2, time1)
00451               << ", write_time: "
00452               << difftime(time3, time2)
00453               << std::endl;
00454   }
00455 }
00456 
00457 EdgeStatus EBMesher::get_edge_status(const double dP, int& iSkip)
00458 {
00459   if (m_nStatus == INSIDE) { // previous inside
00460     if (dP < m_dSecondP) {
00461       m_nMove = 0;
00462       iSkip = m_iSecondP;
00463       return INSIDE;
00464     }
00465     else {
00466       if (is_shared_overlapped_surf(m_iInter - 1)) m_nMove = 1;
00467       else m_nMove = 2;
00468         
00469       // skip shared and overlapped interesections
00470       while (m_vIntersection[m_iInter].distance -
00471              m_vIntersection[m_iInter - 1].distance < 1e-7 &&
00472              (unsigned int) (m_iInter + 1) < m_vIntersection.size()) m_iInter++;        
00473       
00474       return BOUNDARY;
00475     }
00476   }
00477   else if (m_nStatus == OUTSIDE) { // previous outside
00478     if (dP < m_dFirstP) {
00479       m_nMove = 0;
00480       return OUTSIDE;
00481     }
00482     else {
00483       if (dP < m_dSecondP) m_nMove = 0;
00484       else if (is_shared_overlapped_surf(m_iInter - 1)) m_nMove = 1;
00485       else m_nMove = 2;
00486 
00487       // skip shared and overlapped interesections
00488       while (m_vIntersection[m_iInter].distance -
00489              m_vIntersection[m_iInter - 1].distance < 1e-7 &&
00490              (unsigned int) (m_iInter + 1) < m_vIntersection.size()) m_iInter++;
00491 
00492       return BOUNDARY;
00493     }
00494   }
00495   else if (m_nStatus == BOUNDARY) { // previous boundary
00496     if (dP < m_dFirstP) {
00497       m_nMove = 0;
00498       iSkip = m_iFirstP;
00499       return OUTSIDE;
00500     }
00501     else {
00502       if (dP < m_dSecondP) {
00503         m_nMove = 0;
00504         if (m_prevPnt < m_dFirstP) return BOUNDARY;
00505         else {
00506           iSkip = m_iSecondP;
00507           return INSIDE;
00508         }
00509       }
00510       else {
00511         if (is_shared_overlapped_surf(m_iInter - 1)) m_nMove = 1;
00512         else m_nMove = 2;
00513 
00514         // skip shared and overlapped interesections
00515         while (m_vIntersection[m_iInter].distance -
00516                m_vIntersection[m_iInter - 1].distance < 1e-7 &&
00517                (unsigned int) (m_iInter + 1) < m_vIntersection.size()) m_iInter++;
00518 
00519         if (m_prevPnt < m_dSecondP) return BOUNDARY;
00520         else return OUTSIDE;
00521       }
00522     }
00523   }
00524   else {
00525     std::cerr << "Couldn't get edge status." << std::endl;
00526     return INSIDE;
00527   }
00528 }
00529 
00530 bool EBMesher::set_neighbor_hex_status(int dir, int i, int j, int k)
00531 {
00532   int iElem;
00533   int otherDir1 = (dir + 1)%3;
00534   int otherDir2 = (dir + 2)%3;
00535 
00536   if (dir == 0) { // x coordinate ray
00537     m_iElem = j*m_nDiv[0]*m_nDiv[1] + i*m_nDiv[0] + k;
00538     if (!set_hex_status(m_iElem, m_nStatus, dir)) return false;
00539     iElem = m_iElem - m_nDiv[dir];
00540     if (!set_hex_status(iElem, m_nStatus, dir)) return false;
00541     iElem -= m_nDiv[dir]*m_nDiv[otherDir1];
00542     if (!set_hex_status(iElem, m_nStatus, dir)) return false;
00543     iElem += m_nDiv[dir];
00544     if (!set_hex_status(iElem, m_nStatus, dir)) return false;
00545   }
00546   else if (dir == 1) { // y coordinate ray
00547     m_iElem = i*m_nDiv[1]*m_nDiv[0] + k*m_nDiv[0] + j;
00548     if (!set_hex_status(m_iElem, m_nStatus, dir)) return false;
00549     iElem = m_iElem - 1;
00550     if (!set_hex_status(iElem, m_nStatus, dir)) return false;
00551     iElem -= m_nDiv[dir]*m_nDiv[otherDir2];
00552     if (!set_hex_status(iElem++, m_nStatus, dir)) return false;
00553     if (!set_hex_status(iElem, m_nStatus, dir)) return false;
00554   }
00555   else if (dir == 2) { // z coordinate ray
00556     m_iElem = k*m_nDiv[0]*m_nDiv[1] + j*m_nDiv[0] + i;
00557     if (!set_hex_status(m_iElem, m_nStatus, dir)) return false;
00558     iElem = m_iElem - 1;
00559     if (!set_hex_status(iElem, m_nStatus, dir)) return false;
00560     iElem -= m_nDiv[otherDir1];
00561     if (!set_hex_status(iElem++, m_nStatus, dir)) return false;
00562     if (!set_hex_status(iElem, m_nStatus, dir)) return false;
00563   }
00564 
00565   return true;
00566 }
00567 
00568 bool EBMesher::set_hex_status(int index, EdgeStatus value, int dir)
00569 {
00570   if (index < 0 || index > m_nHex - 1) {
00571     return false;
00572   }
00573 
00574   if (m_vnHexStatus[index] != 2) {
00575     if (value == INSIDE) {
00576       m_vnHexStatus[index] = 0;
00577     }
00578     else if (value == OUTSIDE) {
00579       m_vnHexStatus[index] = 1;
00580     }
00581     else if (value == BOUNDARY) {
00582       m_vnHexStatus[index] = 2;
00583     }
00584   }
00585 
00586   return true;
00587 }
00588 
00589 bool EBMesher::set_edge_status(int dir)
00590 {
00591   // set boundary cut information to edge
00592   std::vector<double> vdCutFraction;
00593   if (m_nMove == 0) vdCutFraction.push_back((m_curPnt - m_dFirstP)/m_dIntervalSize[dir] - 1.);
00594   else if (m_nMove == 1) vdCutFraction.push_back(1. - (m_curPnt - m_dSecondP)/m_dIntervalSize[dir]);
00595   else if (m_nMove == 2) {
00596     vdCutFraction.push_back(1. - (m_curPnt - m_dSecondP)/m_dIntervalSize[dir]);
00597     if (m_dFirstP < m_curPnt) {
00598       vdCutFraction.push_back((m_curPnt - m_dFirstP)/m_dIntervalSize[dir] - 1.);
00599     }
00600   }
00601 
00602   std::map<int, CutFraction>::iterator iter = m_mdCutFraction.find(m_iElem);
00603   if (iter == m_mdCutFraction.end()) { // not exist
00604     CutFraction cFraction(dir, vdCutFraction);
00605     m_mdCutFraction[m_iElem] = cFraction;
00606   }
00607   else { // exist
00608     iter->second.add(dir, vdCutFraction);
00609   }
00610 
00611   //m_nFraction += vdCutFraction.size();
00612 
00613   return true;
00614 }
00615 
00616 void EBMesher::set_division()
00617 {
00618   int i, j;
00619   moab::CartVect box_center, box_axis1, box_axis2, box_axis3,
00620     min_cart_box(HUGE_VAL, HUGE_VAL, HUGE_VAL),
00621     max_cart_box(0., 0., 0.);
00622   
00623   moab::ErrorCode rval = m_hObbTree->box(m_hTreeRoot, box_center.array(),
00624                                      box_axis1.array(), box_axis2.array(),
00625                                      box_axis3.array());
00626   MBERRCHK(rval, mk_core()->moab_instance());
00627 
00628   // cartesian box corners
00629   moab::CartVect corners[8] = {box_center + box_axis1 + box_axis2 + box_axis3,
00630                                box_center + box_axis1 + box_axis2 - box_axis3,
00631                                box_center + box_axis1 - box_axis2 + box_axis3,
00632                                box_center + box_axis1 - box_axis2 - box_axis3,
00633                                box_center - box_axis1 + box_axis2 + box_axis3,
00634                                box_center - box_axis1 + box_axis2 - box_axis3,
00635                                box_center - box_axis1 - box_axis2 + box_axis3,
00636                                box_center - box_axis1 - box_axis2 - box_axis3};
00637 
00638   // get the max, min cartesian box corners
00639   for (i = 0; i < 8; i++) {
00640     for (j = 0; j < 3; j++) {
00641       if (corners[i][j] < min_cart_box[j]) min_cart_box[j] = corners[i][j];
00642       if (corners[i][j] > max_cart_box[j]) max_cart_box[j] = corners[i][j];
00643     }
00644   }
00645   moab::CartVect length = max_cart_box - min_cart_box;
00646   
00647   // default value is adjusted to large geometry file
00648   // interval_size_estimate : 2*L*sqrt(2*PI*sqrt(2)/# of tris)
00649   if (m_dInputSize < 0.) {
00650     int n_tri;
00651     rval = moab_instance()->
00652       get_number_entities_by_dimension(reinterpret_cast<moab::EntityHandle>(m_hRootSet),
00653                                        2, n_tri);
00654     MBERRCHK(rval, mk_core()->moab_instance());
00655     
00656     double box_length_ave = (length[0] + length[1] + length[2])/3.;
00657     m_dInputSize = 2.*box_length_ave*sqrt(8.886/n_tri);
00658   }
00659 
00660   for (i = 0; i < 3; i++) {
00661     m_nDiv[i] = length[i]/m_dInputSize;
00662     if (m_nDiv[i] < m_nAddLayer) m_nDiv[i] = m_nAddLayer; // make "m_nAddLayer" elements larger than bounding box
00663     m_dIntervalSize[i] = m_dInputSize;
00664     //if (m_nDiv[i]*.07 > m_nAddLayer) m_nDiv[i] += m_nDiv[i]*.07;
00665     //else m_nDiv[i] += m_nAddLayer;
00666     m_nDiv[i] += m_nAddLayer;
00667     m_nNode[i] = m_nDiv[i] + 1;
00668     m_origin_coords[i] = box_center[i] - .5*m_nDiv[i]*m_dIntervalSize[i];
00669   }
00670 
00671   m_nHex = m_nDiv[0]*m_nDiv[1]*m_nDiv[2];
00672 
00673   std::cout << "# of hex: " << m_nHex << ", interval size: "
00674             << m_dInputSize << std::endl;
00675 
00676   std::cout << "# of division: " << m_nDiv[0] << ","
00677             << m_nDiv[1] << "," << m_nDiv[2] << std::endl;
00678 }
00679 
00680 int EBMesher::set_division(double* min, double* max)
00681 {
00682   for (int i = 0; i < 3; i++) {
00683     m_dIntervalSize[i] = (max[i] - min[i])/m_nDiv[i];
00684     m_nNode[i] = m_nDiv[i] + 1;
00685     m_origin_coords[i] = min[i];
00686   }
00687 
00688   m_nHex = m_nDiv[0]*m_nDiv[1]*m_nDiv[2];
00689 
00690   std::cout << "# of hex: " << m_nHex << ", interval_size: "
00691             << m_dIntervalSize[0] << ", " << m_dIntervalSize[1] << ", "
00692             << m_dIntervalSize[2] << std::endl;
00693 
00694   std::cout << "# of division: " << m_nDiv[0] << ","
00695             << m_nDiv[1] << "," << m_nDiv[2] << std::endl;
00696 
00697   return iBase_SUCCESS;
00698 }
00699 
00700 int EBMesher::make_scd_hexes()
00701 {
00702   // create vertex and hex sequences
00703   int i;
00704   double max_coords[3];
00705   (void) max_coords;
00706   for (i = 0; i < 3; i++) {
00707     max_coords[i] = m_origin_coords[i] + m_nDiv[i]*m_dIntervalSize[i];
00708   }
00709 
00710   moab::HomCoord coord_min(0, 0, 0);
00711   moab::HomCoord coord_max(m_nDiv[0], m_nDiv[1], m_nDiv[2]);
00712   moab::EntitySequence* vertex_seq = NULL;
00713   moab::EntitySequence* cell_seq = NULL;
00714   moab::EntityHandle vs, cs;
00715   moab::Core *mbcore = dynamic_cast<moab::Core*>(moab_instance());
00716 
00717   moab::ErrorCode rval = mbcore->create_scd_sequence(coord_min, coord_max, moab::MBVERTEX, 1, vs, vertex_seq);
00718   MBERRCHK(rval, mk_core()->moab_instance());
00719   
00720   mbcore->create_scd_sequence(coord_min, coord_max, moab::MBHEX, 1, cs, cell_seq);
00721   MBERRCHK(rval, mk_core()->moab_instance());
00722 
00723   moab::HomCoord p2(coord_max.i(), coord_min.j(), coord_min.k());
00724   moab::HomCoord p3(coord_min.i(), coord_max.j(), coord_min.k()); 
00725   
00726   rval = mbcore->add_vsequence(vertex_seq, cell_seq, coord_min, coord_min,
00727                                         p2, p2, p3, p3);
00728   MBERRCHK(rval, mk_core()->moab_instance());
00729 
00730   int nTotNode = m_nNode[0]*m_nNode[1]*m_nNode[2];
00731   int ii, jj, kk;
00732   moab::EntityHandle handle;
00733   double dumv[3];
00734   for (i = 0, handle = vs; i < nTotNode; i++, handle++) {
00735     ii = (i%(m_nNode[0]*m_nNode[1]))%m_nNode[0];
00736     jj = (i%(m_nNode[0]*m_nNode[1]))/m_nNode[0];
00737     kk = i/m_nNode[0]/m_nNode[1];
00738     dumv[0] = m_origin_coords[0] + ii*m_dIntervalSize[0];
00739     dumv[1] = m_origin_coords[1] + jj*m_dIntervalSize[1];
00740     dumv[2] = m_origin_coords[2] + kk*m_dIntervalSize[2];
00741     rval = moab_instance()->set_coords(&handle, 1, dumv);
00742     MBERRCHK(rval, mk_core()->moab_instance());
00743   }
00744 
00745   m_iStartHex = moab_instance()->id_from_handle(cs);
00746 
00747   return iBase_SUCCESS;
00748 }
00749 
00750 iBase_TagHandle EBMesher::get_tag(const char* name, int size,
00751                                      unsigned store, moab::DataType type,
00752                                      const void* def_value,
00753                                      bool create_if_missing) 
00754 {
00755   moab::Tag retval = 0;
00756   /*moab::ErrorCode result = moab_instance()->tag_create(name, size, store, type,
00757                                                    retval, def_value,
00758                                                    create_if_missing);*/
00759   store = store | moab::MB_TAG_CREAT;
00760   if(!create_if_missing)
00761     store = store | moab::MB_TAG_EXCL;
00762   moab::ErrorCode result = moab_instance()->tag_get_handle(name, size, type, retval, store,
00763                                                       def_value);
00764   if (create_if_missing && moab::MB_SUCCESS != result) {
00765     std::cerr << "Couldn't find nor create tag named " << name << std::endl;
00766   }
00767   
00768   return (iBase_TagHandle) retval;
00769 }
00770 
00771 iBase_TagHandle EBMesher::get_various_length_tag(const char* name,
00772                                                     unsigned store, moab::DataType type)
00773 {
00774   moab::Tag retval = 0;
00775   store = store | moab::MB_TAG_VARLEN | moab::MB_TAG_CREAT;
00776   moab::ErrorCode result = moab_instance()->
00777     tag_get_handle( name, 1, type, retval, store);
00778   if (moab::MB_SUCCESS != result) {
00779     std::cerr << "Couldn't find nor create tag named " << name << std::endl;
00780   }
00781   
00782   return (iBase_TagHandle) retval;
00783 }
00784 
00785 void EBMesher::set_tag_info()
00786 {
00787   // get all hexes
00788   int i, j, k;
00789   iBase_ErrorType err;
00790   std::vector<iBase_EntityHandle> hex_handles;
00791   err = mk_core()->imesh_instance()->getEntities(m_hRootSet, iBase_REGION,
00792                                                  iMesh_HEXAHEDRON, hex_handles);
00793   IBERRCHK(err, "Failed to get hexes.\n");
00794 
00795   err = mk_core()->imesh_instance()->setIntArrData(&hex_handles[0], m_nHex,
00796                                                    m_elemStatusTag,
00797                                                    &m_vnHexStatus[0]);
00798   IBERRCHK(err, "Failed to set hex element status data.");
00799 
00800   // set cut fraction info to boundary hexes
00801   int nBndrHex = m_mdCutFraction.size();
00802   std::vector<iBase_EntityHandle> hvBndrHex(nBndrHex);
00803   int* frac_size = new int[nBndrHex];
00804   int* frac_leng = new int[3*nBndrHex];
00805   int nFracSize, nTempSize, iHex, nTolFrac = 0;
00806   int nDoubleSize = sizeof(double);
00807   std::map<int, CutFraction>::iterator iter = m_mdCutFraction.begin();
00808   std::map<int, CutFraction>::iterator end_iter = m_mdCutFraction.end();
00809   for (i = 0; iter != end_iter; iter++, i++) {
00810     iHex = iter->first;
00811     hvBndrHex[i] = hex_handles[iHex];
00812 
00813     nFracSize = 0;
00814     for (j = 0; j < 3; j++) {
00815       nTempSize = iter->second.fraction[j].size();
00816       nFracSize += nTempSize;
00817       frac_leng[3*i + j] = nTempSize;
00818     }
00819     frac_size[i] = nDoubleSize*nFracSize; // sizeof(double)*
00820     nTolFrac += nFracSize;
00821   }
00822 
00823   int iFrac = 0;
00824   double* fractions = new double[nTolFrac];
00825   std::vector<const void*> frac_data_pointer(nBndrHex);
00826   iter = m_mdCutFraction.begin();
00827   for (i = 0; iter != end_iter; iter++, i++) {
00828     for (j = 0; j < 3; j++) {
00829       nFracSize = iter->second.fraction[j].size();
00830       frac_data_pointer[i] = &fractions[iFrac];
00831       for (k = 0; k < nFracSize; k++) {
00832         fractions[iFrac++] = iter->second.fraction[j][k];
00833       }
00834     }
00835   }
00836 
00837   err = mk_core()->imesh_instance()->setIntArrData(&hvBndrHex[0], nBndrHex,
00838                                                    m_edgeCutFracLengthTag,
00839                                                    frac_leng);
00840   IBERRCHK(err, "Failed to set cut fraction sizes to hex.");
00841 
00842   moab::ErrorCode rval = moab_instance()->tag_set_by_ptr(reinterpret_cast<moab::Tag> (m_edgeCutFracTag),
00843                                                    reinterpret_cast<moab::EntityHandle*> (&hvBndrHex[0]),
00844                                                    nBndrHex, &frac_data_pointer[0], frac_size);
00845   MBERRCHK(rval, mk_core()->moab_instance());
00846   
00847   delete [] frac_size;
00848   delete [] frac_leng;
00849   delete [] fractions;
00850 }
00851 
00852 void EBMesher::fire_rays(int dir)
00853 {
00854   // ray fire
00855   int i, j, k, l, index[3];
00856   double tolerance = 1e-12;
00857   double rayLength = m_nDiv[dir]*m_dIntervalSize[dir];
00858   int iNodeStart, iNodeEnd, nIntersect, nNodeSlice;
00859   double startPnt[3], endPnt[3];
00860   int otherDir1 = (dir + 1)%3;
00861   int otherDir2 = (dir + 2)%3;
00862  // harmless as this expression (does nothing) so that a compiler sees it is used.
00863   (void) iNodeEnd;
00864   (void) nNodeSlice;
00865 
00866   for (j = 1; j < m_nNode[otherDir2] - 1; j++) {
00867     for (i = 1; i < m_nNode[otherDir1] - 1; i++) {
00868 
00869       // get ray start and end points
00870       if (dir == 0) { // x coordinate ray
00871         iNodeStart = j*m_nNode[dir]*m_nNode[otherDir1] + i*m_nNode[dir];
00872         iNodeEnd = iNodeStart + m_nNode[dir] - 1;
00873         nNodeSlice = 1;
00874         index[0] = 0;
00875         index[1] = i;
00876         index[2] = j;
00877       }
00878       else if (dir == 1) { // y coordinate ray
00879         iNodeStart = i*m_nNode[otherDir2]*m_nNode[dir] + j;
00880         iNodeEnd = iNodeStart + m_nNode[otherDir2]*(m_nNode[dir] - 1);
00881         nNodeSlice = m_nNode[otherDir2];
00882         index[0] = j;
00883         index[1] = 0;
00884         index[2] = i;
00885       }
00886       else if (dir == 2) { // z coordinate ray
00887         iNodeStart = j*m_nNode[otherDir1] + i;
00888         iNodeEnd = iNodeStart + m_nNode[otherDir1]*m_nNode[otherDir2]*(m_nNode[dir] - 1);
00889         nNodeSlice = m_nNode[otherDir1]*m_nNode[otherDir2];
00890         index[0] = i;
00891         index[1] = j;
00892         index[2] = 0;
00893       }
00894       else IBERRCHK(iBase_FAILURE, "Ray direction should be 0 to 2.");
00895       
00896       for (l = 0; l < 3; l++) {
00897         if (l == dir) {
00898           startPnt[l] = m_origin_coords[l];
00899           endPnt[l] = m_origin_coords[l] + m_nDiv[dir]*m_dIntervalSize[l];
00900         }
00901         else {
00902           startPnt[l] = m_origin_coords[l] + index[l]*m_dIntervalSize[l];
00903           endPnt[l] = startPnt[l];
00904         }
00905       }
00906       
00907       k = 0;
00908       m_iInter = 0;
00909       if (!fire_ray(nIntersect, startPnt, endPnt, tolerance, dir, rayLength)) {
00910         IBERRCHK(iBase_FAILURE, "Couldn't fire ray.");
00911       }
00912       
00913       if (nIntersect > 0) {
00914         m_iFirstP = m_vIntersection[m_iInter].distance/m_dIntervalSize[dir];
00915         m_dFirstP = startPnt[dir] + m_vIntersection[m_iInter++].distance;
00916         m_iSecondP = m_vIntersection[m_iInter].distance/m_dIntervalSize[dir];
00917         m_dSecondP = startPnt[dir] + m_vIntersection[m_iInter++].distance;
00918 
00919         // set outside before the first hit
00920         for (k = 0; k < m_iFirstP; k++) {
00921           m_nStatus = OUTSIDE;
00922           
00923           if (!set_neighbor_hex_status(dir, i, j, k)) {
00924             IBERRCHK(iBase_FAILURE, "Couldn't set neighbor hex status.");
00925           }
00926         }
00927         
00928         for (; k < m_nNode[dir] - 1;) {
00929           int i_skip = 0;
00930           m_curPnt = startPnt[dir] + (k + 1)*m_dIntervalSize[dir];
00931           EdgeStatus preStat = m_nStatus;
00932           m_nStatus = get_edge_status(m_curPnt, i_skip);
00933           m_vnEdgeStatus[dir][i*m_nDiv[dir] + j*m_nDiv[dir]*m_nDiv[otherDir1] + k] = m_nStatus;
00934 
00935           // set status of all hexes sharing the edge
00936           if (!set_neighbor_hex_status(dir, i, j, k)) {
00937             IBERRCHK(iBase_FAILURE, "Couldn't set neighbor hex status.");
00938           }
00939 
00940           if (m_nMove > 0) {
00941             if (m_iInter < nIntersect) {
00942               if (!move_intersections(dir, nIntersect, startPnt)) {
00943                 IBERRCHK(iBase_FAILURE, "Couldn't move intersections.");
00944               }
00945             }
00946             else {
00947               m_nMove = 1;
00948               if (m_nStatus == BOUNDARY && !set_edge_status(dir)) {
00949                 IBERRCHK(iBase_FAILURE, "Couldn't set edge status.");
00950               }
00951               k++;
00952               break; // rest is all outside
00953             }
00954           }
00955           else if (m_nStatus == BOUNDARY && !set_edge_status(dir)) {
00956             IBERRCHK(iBase_FAILURE, "Couldn't set edge status.");
00957           }
00958           else if (m_nStatus == OUTSIDE && preStat == BOUNDARY) { // set outside until next hit
00959             k++;
00960             for (; k < i_skip; k++) {
00961               m_nStatus = OUTSIDE;
00962               if (!set_neighbor_hex_status(dir, i, j, k)) {
00963                 IBERRCHK(iBase_FAILURE, "Couldn't set neighbor hex status.");
00964               }
00965             }
00966           }
00967 
00968           // set cut-cell edge status
00969           if (i_skip > 0) {
00970             m_prevPnt = startPnt[dir] + i_skip*m_dIntervalSize[dir];
00971             k = i_skip;
00972           }
00973           else {
00974             m_prevPnt = m_curPnt;
00975             k++;
00976           }
00977         }
00978       }
00979         
00980       // the rest are all outside
00981       for (; k < m_nNode[dir] - 1; k++) {
00982         m_nStatus = OUTSIDE;
00983         if (!set_neighbor_hex_status(dir, i, j, k)) {
00984           IBERRCHK(iBase_FAILURE, "Couldn't set neighbor hex status.");
00985         }
00986       }
00987     }
00988   }
00989 }
00990 
00991 bool EBMesher::fire_ray(int& nIntersect, double startPnt[3],
00992                       double endPnt[3], double tol, int dir,
00993                       double rayLength)
00994 {
00995   m_vIntersection.clear();
00996   m_vhInterSurf.clear();
00997   m_vhInterFacet.clear();
00998   m_mhOverlappedSurf.clear();
00999   std::vector<double> temp_intersects;
01000   moab::ErrorCode rVal;
01001   moab::OrientedBoxTreeTool::TrvStats stats;
01002   if (m_bUseGeom) { // geometry input
01003     rVal = m_hObbTree->ray_intersect_sets(temp_intersects, m_vhInterSurf,
01004                                           m_vhInterFacet, m_hTreeRoot, tol,
01005                                           startPnt, rayDir[dir], &rayLength,
01006                                           &stats);
01007   }
01008   else { // facet input
01009     std::vector<moab::EntityHandle> dum_facets_out;
01010     rVal = m_hObbTree->ray_intersect_triangles(temp_intersects, dum_facets_out,
01011                                                m_hTreeRoot, tol,
01012                                                startPnt, rayDir[dir], &rayLength,
01013                                                &stats);
01014   }
01015   
01016   nIntersect = temp_intersects.size();
01017   if (moab::MB_SUCCESS != rVal) {
01018     std::cerr << "Failed : ray-triangle intersection." << std::endl;
01019     return false;
01020   }
01021   
01022   if (nIntersect > 0) {
01023     m_vIntersection.resize(nIntersect);
01024     
01025     for (int l = 0; l < nIntersect; l++) {
01026       IntersectDist temp_inter_dist(temp_intersects[l], l);
01027       m_vIntersection[l] = temp_inter_dist;
01028     }
01029     std::sort(m_vIntersection.begin(), m_vIntersection.end(), less_intersect);
01030 
01031     if (nIntersect > 1) {
01032       bool bMoveOnce;
01033       m_nIteration = 0;
01034       m_iOverlap = 0;
01035       
01036       if (m_bUseGeom) { // if there are geometry info
01037         remove_intersection_duplicates();
01038       }
01039       else if (is_ray_move_and_set_overlap_surf(bMoveOnce)) { // facet geom case
01040         if (!move_ray(nIntersect, startPnt, endPnt, tol, dir, bMoveOnce)) {
01041           std::cerr << "Number of Intersection between edges and ray should be even." << std::endl;
01042           return false;
01043         }
01044       }
01045       nIntersect = m_vIntersection.size();
01046     }
01047   }
01048   
01049   return true;
01050 }
01051 
01052 bool EBMesher::move_intersections(int n_dir, int n_inter, double start_pnt[3])
01053 {
01054   if (m_nMove > 0) {
01055     if (m_iInter < n_inter) {
01056       if (m_nMove == 1) {
01057         do {
01058           if (m_nStatus == BOUNDARY && !set_edge_status(n_dir)) return false;
01059           m_iFirstP = m_iSecondP;
01060           m_dFirstP = m_dSecondP;
01061           m_iSecondP = m_vIntersection[m_iInter].distance/m_dIntervalSize[n_dir];
01062           m_dSecondP = start_pnt[n_dir] + m_vIntersection[m_iInter++].distance;
01063         }
01064         while (m_dSecondP < m_curPnt && m_iInter < n_inter);
01065       }
01066       else if (m_nMove == 2) {
01067         do {
01068           m_iFirstP = m_vIntersection[m_iInter].distance/m_dIntervalSize[n_dir];
01069           m_dFirstP = start_pnt[n_dir] + m_vIntersection[m_iInter++].distance;
01070           if (m_nStatus == BOUNDARY && !set_edge_status(n_dir)) return false;
01071           if (m_iInter < n_inter) {
01072             m_iSecondP = m_vIntersection[m_iInter].distance/m_dIntervalSize[n_dir];
01073             m_dSecondP = start_pnt[n_dir] + m_vIntersection[m_iInter++].distance;
01074           }
01075           else {
01076             m_iSecondP = m_iFirstP;
01077             m_dSecondP = m_dFirstP;
01078           }
01079         }
01080         while (m_dSecondP < m_curPnt && m_iInter < n_inter);
01081       }
01082     }
01083   }
01084 
01085   return true;
01086 }
01087 
01088 bool EBMesher::is_shared_overlapped_surf(int index)
01089 {
01090   int nParent;
01091   iBase_ErrorType err;
01092   moab::EntityHandle hSurf;
01093   if (m_bUseGeom) {
01094     hSurf = m_vhInterSurf[m_vIntersection[index].index];
01095     err = mk_core()->imesh_instance()->getNumPrnt(reinterpret_cast<iBase_EntitySetHandle> (hSurf),
01096                                                   1, nParent);
01097     IBERRCHK(err, "Failed to get number of surface parents.\n");
01098 
01099     if (nParent > 1) return true;
01100   }
01101   else hSurf = m_vIntersection[index].index;
01102 
01103   return m_mhOverlappedSurf.count(hSurf) > 0;
01104 }
01105 
01106 void EBMesher::get_grid_and_edges_techX(double* boxMin, double* boxMax, int* nDiv,
01107                                       std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellSurfEdge,
01108                                       std::vector<int>& rvnInsideCell, bool isCornerExterior)
01109 {
01110   int i, j, err, ii, jj, kk, iHex;
01111   for (i = 0; i < 3; i++) {
01112     boxMin[i] = m_origin_coords[i];
01113     boxMax[i] = m_origin_coords[i] + m_dIntervalSize[i]*m_nDiv[i];
01114     nDiv[i] = m_nDiv[i];
01115   }
01116   
01117   // get all hexes
01118   std::vector<iBase_EntityHandle> hex_handles;
01119   err = mk_core()->imesh_instance()->getEntities(m_hRootSet, iBase_REGION,
01120                                                  iMesh_HEXAHEDRON, hex_handles);
01121   IBERRCHK(err, "Failed to get hexes.\n");
01122   
01123   // get hex status
01124   int hex_size = hex_handles.size();
01125   /*
01126   //int* hex_status = new int[hex_size];
01127   //int* hex_status;
01128   //std::vector<int> hex_status(hex_size);
01129   int temp = 0;
01130   int* hex_status = &temp;
01131   err = mk_core()->imesh_instance()->getIntArrData(&hex_handles[0], hex_size,
01132                                                    m_elemStatusTag, hex_status);
01133   */
01134   int error;
01135   int* hex_status = new int[hex_size];
01136   int hex_status_alloc = sizeof(int)*hex_size;
01137   int hex_status_size = 0;
01138   iMesh_getIntArrData(m_mesh, &hex_handles[0],
01139                       hex_size, m_elemStatusTag,
01140                       &hex_status, &hex_status_alloc,
01141                       &hex_status_size, &error);
01142   IBERRCHK(error, "Failed to get hex status.");
01143 
01144   // get inside and boundary hexes
01145   int nInside = 0;
01146   int nOutside = 0;
01147   int nBoundary = 0;
01148   rvnInsideCell.clear();
01149   for (i = 0; i < hex_size; i++) {
01150     if (hex_status[i] == 0) { // if inside
01151       iHex = moab_instance()->id_from_handle(reinterpret_cast<moab::EntityHandle>
01152                                              (hex_handles[i])) - m_iStartHex;
01153       rvnInsideCell.push_back((iHex%(m_nDiv[0]*m_nDiv[1]))%m_nDiv[0]);
01154       rvnInsideCell.push_back((iHex%(m_nDiv[0]*m_nDiv[1]))/m_nDiv[0]);
01155       rvnInsideCell.push_back(iHex/m_nDiv[0]/m_nDiv[1]);
01156       nInside++;
01157     }
01158     else if (hex_status[i] == 1) nOutside++;
01159     else if (hex_status[i] == 2) nBoundary++;
01160     else IBERRCHK(err, "Element status should be one of inside/outside/boundary."); 
01161   }
01162   std::cout << "# of inside, outside, boundary elem: " << nInside
01163             << ", " << nOutside << ", " << nBoundary << std::endl;
01164   delete [] hex_status;
01165 
01166   // get cut-cell fractions
01167   double dFrac[4];
01168   std::map<int, CutFraction>::iterator iter = m_mdCutFraction.begin();
01169   std::map<int, CutFraction>::iterator end_iter = m_mdCutFraction.end();
01170   
01171   for (; iter != end_iter; iter++) { // for each cut-cell
01172     // get i, j, k index from handle
01173     iHex = iter->first;
01174     ii = (iHex%(m_nDiv[0]*m_nDiv[1]))%m_nDiv[0];
01175     jj = (iHex%(m_nDiv[0]*m_nDiv[1]))/m_nDiv[0];
01176     kk = iHex/m_nDiv[0]/m_nDiv[1];
01177 
01178     // surface 1
01179     CutFraction cutFrac = iter->second;
01180     if (cutFrac.fraction[1].size() > 0) {
01181       dFrac[0] = cutFrac.fraction[1][0];
01182       if (dFrac[0] < 0.) dFrac[0] *= -1.;
01183     }
01184     else dFrac[0] = -1.;
01185     if (cutFrac.fraction[2].size() > 0) {
01186       dFrac[3] = cutFrac.fraction[2][0];
01187       if (dFrac[3] < 0.) dFrac[3] *= -1.;
01188     }
01189     else dFrac[3] = -1.;
01190     dFrac[1] = get_edge_fraction(iHex + m_nDiv[0], 2);
01191     dFrac[2] = get_edge_fraction(iHex + m_nDiv[0]*m_nDiv[1], 1);
01192     
01193     if (dFrac[0] > 0. || dFrac[1] > 0. || dFrac[2] > 0. || dFrac[3] > 0.) { // if surface is cut
01194       CutCellSurfEdgeKey surfKey(ii, jj, kk, 0);
01195       std::vector<double> edges;
01196       if (dFrac[0] < 0.) dFrac[0] = get_uncut_edge_fraction(ii, jj, kk, 1);
01197       if (dFrac[1] < 0.) dFrac[1] = get_uncut_edge_fraction(ii, jj + 1, kk, 2);
01198       if (dFrac[2] < 0.) dFrac[2] = get_uncut_edge_fraction(ii, jj, kk + 1, 1);
01199       if (dFrac[3] < 0.) dFrac[3] = get_uncut_edge_fraction(ii, jj, kk, 2);
01200       for (j = 0; j < 4; j++) edges.push_back(dFrac[j]);
01201       rmdCutCellSurfEdge[surfKey] = edges;
01202     }
01203 
01204     // surface 2
01205     cutFrac = iter->second;
01206     if (cutFrac.fraction[0].size() > 0) {
01207       dFrac[0] = cutFrac.fraction[0][0];
01208       if (dFrac[0] < 0.) dFrac[0] *= -1.;
01209     }
01210     else dFrac[0] = -1.;
01211     if (cutFrac.fraction[1].size() > 0) {
01212       dFrac[3] = cutFrac.fraction[1][0];
01213       if (dFrac[3] < 0.) dFrac[3] *= -1.;
01214     }
01215     else dFrac[3] = -1.;
01216     dFrac[1] = get_edge_fraction(iHex + 1, 2);
01217     dFrac[2] = get_edge_fraction(iHex + m_nDiv[0]*m_nDiv[1], 0);
01218     
01219     if (dFrac[0] > 0. || dFrac[1] > 0. || dFrac[2] > 0. || dFrac[3] > 0.) { // if surface is cut
01220       CutCellSurfEdgeKey surfKey(ii, jj, kk, 1);
01221       std::vector<double> edges;
01222       if (dFrac[0] < 0.) dFrac[0] = get_uncut_edge_fraction(ii, jj, kk, 0);
01223       if (dFrac[1] < 0.) dFrac[1] = get_uncut_edge_fraction(ii + 1, jj, kk, 2);
01224       if (dFrac[2] < 0.) dFrac[2] = get_uncut_edge_fraction(ii, jj, kk + 1, 0);
01225       if (dFrac[3] < 0.) dFrac[3] = get_uncut_edge_fraction(ii, jj, kk, 2);
01226       for (j = 0; j < 4; j++) edges.push_back(dFrac[j]);
01227       rmdCutCellSurfEdge[surfKey] = edges;
01228     }
01229     
01230     // surface 3
01231     cutFrac = iter->second;
01232     if (cutFrac.fraction[0].size() > 0) {
01233       dFrac[0] = cutFrac.fraction[0][0];
01234       if (dFrac[0] < 0.) dFrac[0] *= -1.;
01235     }
01236     else dFrac[0] = -1.;
01237     if (cutFrac.fraction[1].size() > 0) {
01238       dFrac[3] = cutFrac.fraction[1][0];
01239       if (dFrac[3] < 0.) dFrac[3] *= -1.;
01240     }
01241     else dFrac[3] = -1.;
01242     dFrac[1] = get_edge_fraction(iHex + 1, 1);
01243     dFrac[2] = get_edge_fraction(iHex + m_nDiv[0], 0);
01244 
01245     if (dFrac[0] > 0. || dFrac[1] > 0. || dFrac[2] > 0. || dFrac[3] > 0.) { // if surface is cut
01246       CutCellSurfEdgeKey surfKey(ii, jj, kk, 2);
01247       std::vector<double> edges;
01248       if (dFrac[0] < 0.) dFrac[0] = get_uncut_edge_fraction(ii, jj, kk, 0);
01249       if (dFrac[1] < 0.) dFrac[1] = get_uncut_edge_fraction(ii + 1, jj, kk, 1);
01250       if (dFrac[2] < 0.) dFrac[2] = get_uncut_edge_fraction(ii, jj + 1, kk, 0);
01251       if (dFrac[3] < 0.) dFrac[3] = get_uncut_edge_fraction(ii, jj, kk, 1);
01252       for (j = 0; j < 4; j++) edges.push_back(dFrac[j]);
01253       rmdCutCellSurfEdge[surfKey] = edges;
01254     }
01255   }
01256   
01257   if (debug_ebmesh) export_fraction_edges(rmdCutCellSurfEdge);
01258 }
01259 
01260 bool EBMesher::get_grid_and_edges(double* boxMin, double* boxMax, int* nDiv,
01261                                      std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellEdge,
01262                                      std::vector<int>& rvnInsideCell, bool isCornerExterior)
01263 {
01264   int i, ii, jj, kk, iHex;
01265   for (i = 0; i < 3; i++) {
01266     boxMin[i] = m_origin_coords[i];
01267     boxMax[i] = m_origin_coords[i] + m_dIntervalSize[i]*m_nDiv[i];
01268     nDiv[i] = m_nDiv[i];
01269   }
01270   
01271   if (!get_inside_hex(rvnInsideCell)) return false;
01272 
01273   // get cut-cell fractions
01274   std::map<int, CutFraction>::iterator iter = m_mdCutFraction.begin();
01275   std::map<int, CutFraction>::iterator end_iter = m_mdCutFraction.end();
01276   
01277   for (; iter != end_iter; iter++) { // for each cut-cell
01278     // get i, j, k index from handle
01279     iHex = iter->first;
01280     ii = (iHex%(m_nDiv[0]*m_nDiv[1]))%m_nDiv[0];
01281     jj = (iHex%(m_nDiv[0]*m_nDiv[1]))/m_nDiv[0];
01282     kk = iHex/m_nDiv[0]/m_nDiv[1];
01283 
01284     for (i = 0; i < 3; i++) {
01285       std::vector<double> fractions(iter->second.fraction[i]);
01286       CutCellSurfEdgeKey edgeKey(ii, jj, kk, i);
01287       rmdCutCellEdge[edgeKey] = fractions;
01288     }
01289   }
01290 
01291   if (debug_ebmesh) return export_fraction_points(rmdCutCellEdge);
01292  
01293   return true;
01294 }
01295 
01296 bool EBMesher::get_inside_hex(std::vector<int>& rvnInsideCell)
01297 {
01298   int i, err, iHex;
01299   
01300   // get all hexes
01301   std::vector<iBase_EntityHandle> hex_handles;
01302   err = mk_core()->imesh_instance()->getEntities(m_hRootSet, iBase_REGION,
01303                                                  iMesh_HEXAHEDRON, hex_handles);
01304   IBERRCHK(err, "Failed to get hexes.\n");
01305   
01306   // get hex status
01307   int hex_size = hex_handles.size();
01308   /*
01309   //int* hex_status = new int[hex_size];
01310   //int* hex_status;
01311   std::vector<int> hex_status(hex_size);
01312   err = mk_core()->imesh_instance()->getIntArrData(&hex_handles[0], hex_size,
01313                                                    m_elemStatusTag, &hex_status[0]);
01314   */
01315   int error;
01316   int* hex_status = new int[hex_size];
01317   int hex_status_alloc = sizeof(int)*hex_size;
01318   int hex_status_size = 0;
01319   iMesh_getIntArrData(m_mesh, &hex_handles[0],
01320                       hex_size, m_elemStatusTag,
01321                       &hex_status, &hex_status_alloc,
01322                       &hex_status_size, &error);
01323   ERRORRF("Failed to get hex status.\n");
01324 
01325   // get inside and boundary hexes
01326   int nInside = 0;
01327   int nOutside = 0;
01328   int nBoundary = 0;
01329   rvnInsideCell.clear();
01330   for (i = 0; i < hex_size; i++) {
01331     if (hex_status[i] == 0) { // if inside
01332       iHex = moab_instance()->id_from_handle(reinterpret_cast<moab::EntityHandle>
01333                                              (hex_handles[i])) - m_iStartHex;
01334       rvnInsideCell.push_back((iHex%(m_nDiv[0]*m_nDiv[1]))%m_nDiv[0]);
01335       rvnInsideCell.push_back((iHex%(m_nDiv[0]*m_nDiv[1]))/m_nDiv[0]);
01336       rvnInsideCell.push_back(iHex/m_nDiv[0]/m_nDiv[1]);
01337       nInside++;
01338     }
01339     else if (hex_status[i] == 1) nOutside++;
01340     else if (hex_status[i] == 2) nBoundary++;
01341     else ERRORRF("Element status should be one of inside/outside/boundary.\n"); 
01342   }
01343   std::cout << "# of inside, outside, boundary elem: " << nInside
01344             << ", " << nOutside << ", " << nBoundary << std::endl;
01345 
01346   delete [] hex_status;
01347   return true;
01348 }
01349 
01350 bool EBMesher::export_fraction_edges(std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& mdCutCellSurfEdge)
01351 {
01352   // export fractions as edge
01353   double curPnt[3], ePnt[6];
01354   std::vector<iBase_EntityHandle> edge_handles;
01355   std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >::iterator itr = mdCutCellSurfEdge.begin();
01356   std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >::iterator e_itr = mdCutCellSurfEdge.end();
01357   for (; itr != e_itr; itr++) {
01358     curPnt[0] = m_origin_coords[0] + itr->first.i*m_dIntervalSize[0];
01359     curPnt[1] = m_origin_coords[1] + itr->first.j*m_dIntervalSize[1];
01360     curPnt[2] = m_origin_coords[2] + itr->first.k*m_dIntervalSize[2];
01361     std::vector<double> edges = itr->second;
01362     
01363     if (itr->first.l == 0) {
01364       if (edges[0] > 0. && edges[0] < 1.) {
01365         ePnt[0] = curPnt[0];
01366         ePnt[1] = curPnt[1];
01367         ePnt[2] = curPnt[2];
01368         ePnt[3] = ePnt[0];
01369         ePnt[4] = ePnt[1] + edges[0]*m_dIntervalSize[1];
01370         ePnt[5] = ePnt[2];
01371         if (!make_edge(ePnt, edge_handles)) return false;
01372       }
01373       if (edges[1] > 0. && edges[1] < 1.) {
01374         ePnt[0] = curPnt[0];
01375         ePnt[1] = curPnt[1] + m_dIntervalSize[1];
01376         ePnt[2] = curPnt[2];
01377         ePnt[3] = ePnt[0];
01378         ePnt[4] = ePnt[1];
01379         ePnt[5] = ePnt[2] + edges[1]*m_dIntervalSize[2];
01380         if (!make_edge(ePnt, edge_handles)) return false;
01381       }
01382       if (edges[2] > 0. && edges[2] < 1.) {
01383         ePnt[0] = curPnt[0];
01384         ePnt[1] = curPnt[1];
01385         ePnt[2] = curPnt[2] + m_dIntervalSize[2];
01386         ePnt[3] = ePnt[0];
01387         ePnt[4] = ePnt[1] + edges[2]*m_dIntervalSize[1];
01388         ePnt[5] = ePnt[2];
01389         if (!make_edge(ePnt, edge_handles)) return false;
01390       }
01391       if (edges[3] > 0. && edges[3] < 1.) {
01392         ePnt[0] = curPnt[0];
01393         ePnt[1] = curPnt[1];
01394         ePnt[2] = curPnt[2];
01395         ePnt[3] = ePnt[0];
01396         ePnt[4] = ePnt[1];
01397         ePnt[5] = ePnt[2] + edges[3]*m_dIntervalSize[2];
01398         if (!make_edge(ePnt, edge_handles)) return false;
01399       }
01400     }
01401     if (itr->first.l == 1) {
01402       if (edges[0] > 0. && edges[0] < 1.) {
01403         ePnt[0] = curPnt[0];
01404         ePnt[1] = curPnt[1];
01405         ePnt[2] = curPnt[2];
01406         ePnt[3] = ePnt[0] + edges[0]*m_dIntervalSize[0];
01407         ePnt[4] = ePnt[1];
01408         ePnt[5] = ePnt[2];
01409         if (!make_edge(ePnt, edge_handles)) return false;
01410       }
01411       if (edges[1] > 0. && edges[1] < 1.) {
01412         ePnt[0] = curPnt[0] + m_dIntervalSize[0];
01413         ePnt[1] = curPnt[1];
01414         ePnt[2] = curPnt[2];
01415         ePnt[3] = ePnt[0];
01416         ePnt[4] = ePnt[1];
01417         ePnt[5] = ePnt[2] + edges[1]*m_dIntervalSize[2];
01418         if (!make_edge(ePnt, edge_handles)) return false;
01419       }
01420       if (edges[2] > 0. && edges[2] < 1.) {
01421         ePnt[0] = curPnt[0];
01422         ePnt[1] = curPnt[1];
01423         ePnt[2] = curPnt[2] + m_dIntervalSize[2];
01424         ePnt[3] = ePnt[0] + edges[2]*m_dIntervalSize[0];
01425         ePnt[4] = ePnt[1];
01426         ePnt[5] = ePnt[2];
01427         if (!make_edge(ePnt, edge_handles)) return false;
01428       }
01429       if (edges[3] > 0. && edges[3] < 1.) {
01430         ePnt[0] = curPnt[0];
01431         ePnt[1] = curPnt[1];
01432         ePnt[2] = curPnt[2];
01433         ePnt[3] = ePnt[0];
01434         ePnt[4] = ePnt[1];
01435         ePnt[5] = ePnt[2] + edges[3]*m_dIntervalSize[2];
01436         if (!make_edge(ePnt, edge_handles)) return false;
01437       }
01438     }
01439     if (itr->first.l == 2) {
01440       if (edges[0] > 0. && edges[0] < 1.) {
01441         ePnt[0] = curPnt[0];
01442         ePnt[1] = curPnt[1];
01443         ePnt[2] = curPnt[2];
01444         ePnt[3] = ePnt[0] + edges[0]*m_dIntervalSize[0];
01445         ePnt[4] = ePnt[1];
01446         ePnt[5] = ePnt[2];
01447         if (!make_edge(ePnt, edge_handles)) return false;
01448       }
01449       if (edges[1] > 0. && edges[1] < 1.) {
01450         ePnt[0] = curPnt[0] + m_dIntervalSize[0];
01451         ePnt[1] = curPnt[1];
01452         ePnt[2] = curPnt[2];
01453         ePnt[3] = ePnt[0];
01454         ePnt[4] = ePnt[1] + edges[1]*m_dIntervalSize[1];
01455         ePnt[5] = ePnt[2];
01456         if (!make_edge(ePnt, edge_handles)) return false;
01457       }
01458       if (edges[2] > 0. && edges[2] < 1.) {
01459         ePnt[0] = curPnt[0];
01460         ePnt[1] = curPnt[1] + m_dIntervalSize[2];
01461         ePnt[2] = curPnt[2];
01462         ePnt[3] = ePnt[0] + edges[2]*m_dIntervalSize[0];
01463         ePnt[4] = ePnt[1];
01464         ePnt[5] = ePnt[2];
01465         if (!make_edge(ePnt, edge_handles)) return false;
01466       }
01467       if (edges[3] > 0. && edges[3] < 1.) {
01468         ePnt[0] = curPnt[0];
01469         ePnt[1] = curPnt[1];
01470         ePnt[2] = curPnt[2];
01471         ePnt[3] = ePnt[0];
01472         ePnt[4] = ePnt[1] + edges[3]*m_dIntervalSize[1];
01473         ePnt[5] = ePnt[2];
01474         if (!make_edge(ePnt, edge_handles)) return false;
01475       }
01476     }
01477   }
01478 
01479   int is_list = 1;
01480   iBase_EntitySetHandle set;
01481   iBase_ErrorType err;
01482   err = mk_core()->imesh_instance()->createEntSet(is_list, set);
01483   IBERRCHK(err, "Couldn't create set.");
01484   
01485   err = mk_core()->imesh_instance()->addEntArrToSet(&edge_handles[0],
01486                                                     edge_handles.size(), set);
01487   IBERRCHK(err, "Couldn't add edges to set.");
01488   
01489   moab::ErrorCode rval = moab_instance()->write_mesh("edges.vtk",
01490                                                      (const moab::EntityHandle*) &set, 1);
01491   MBERRCHK(rval, mk_core()->moab_instance());
01492 
01493   return true;
01494 }
01495 
01496 bool EBMesher::export_fraction_points(std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& mdCutCellEdge)
01497 {
01498   // export fractions as edge
01499   int i, j, dir, nFrc;
01500   iBase_ErrorType err;
01501   double curPnt[3], fracPnt[3], frac;
01502   iBase_EntityHandle v_handle;
01503   std::vector<iBase_EntityHandle> vertex_handles;
01504   std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >::iterator itr = mdCutCellEdge.begin();
01505   std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >::iterator e_itr = mdCutCellEdge.end();
01506   for (; itr != e_itr; itr++) {
01507     curPnt[0] = m_origin_coords[0] + itr->first.i*m_dIntervalSize[0];
01508     curPnt[1] = m_origin_coords[1] + itr->first.j*m_dIntervalSize[1];
01509     curPnt[2] = m_origin_coords[2] + itr->first.k*m_dIntervalSize[2];
01510     dir = itr->first.l;
01511     nFrc = itr->second.size();
01512     
01513     for (i = 0; i < nFrc; i++) {
01514       frac = itr->second[i];
01515       if (frac > 0. && frac < 1.) {
01516         for (j = 0; j < 3; j++) {
01517           if (j == dir) fracPnt[j] = curPnt[j] + frac*m_dIntervalSize[j];
01518           else fracPnt[j] = curPnt[j];
01519         }
01520         err = mk_core()->imesh_instance()->createVtx(fracPnt[0], fracPnt[1],
01521                                                      fracPnt[2], v_handle);
01522         IBERRCHK(err, "Couldn't create vertex.");
01523         vertex_handles.push_back(v_handle);
01524       }
01525     }
01526   }
01527     
01528   int is_list = 1;
01529   moab::ErrorCode result;
01530   iBase_EntitySetHandle set;
01531   err = mk_core()->imesh_instance()->createEntSet(is_list, set);
01532   IBERRCHK(err, "Couldn't create set.");
01533   
01534   err = mk_core()->imesh_instance()->addEntArrToSet(&vertex_handles[0],
01535                                                     vertex_handles.size(), set);
01536   IBERRCHK(err, "Couldn't add vertices to set.");
01537   
01538   result = moab_instance()->write_mesh("frac_vertices.vtk",
01539                                        (const moab::EntityHandle*) &set, 1);
01540   if (moab::MB_SUCCESS != result) {
01541     std::cerr << "Failed to write fraction vertices." << std::endl;
01542     return false;
01543   }
01544 
01545   return true;
01546 }
01547 
01548 bool EBMesher::make_edge(double ePnt[6], std::vector<iBase_EntityHandle>& edge_handles)
01549 {
01550   iBase_ErrorType err;
01551   iBase_EntityHandle vertex_handle[2], edge_handle;
01552   iBase_EntityHandle* pVertexHandle = &vertex_handle[0];
01553   err = mk_core()->imesh_instance()->createVtxArr(2, iBase_INTERLEAVED,
01554                                                   ePnt, pVertexHandle);
01555   IBERRCHK(err, "Failed to create vertices.\n");
01556 
01557   err = mk_core()->imesh_instance()->createEnt(iMesh_LINE_SEGMENT,
01558                                                &vertex_handle[0], 2,
01559                                                edge_handle);
01560   IBERRCHK(err, "Failed to create edge.\n");
01561 
01562   edge_handles.push_back(edge_handle);
01563 
01564   return true;
01565 }
01566 
01567 double EBMesher::get_edge_fraction(int idHex, int dir)
01568 {
01569   std::map<int, CutFraction>::iterator end_iter = m_mdCutFraction.end();
01570   std::map<int, CutFraction>::iterator iter = m_mdCutFraction.find(idHex);
01571   if (iter != end_iter && iter->second.fraction[dir].size() > 0) {
01572     double frac = iter->second.fraction[dir][0];
01573     if (frac < 0.) frac *= -1.;
01574     return frac;
01575   }
01576   else return -1.;
01577 }
01578 
01579 double EBMesher::get_uncut_edge_fraction(int i, int j, int k, int dir)
01580 {
01581   int iEdge;
01582 
01583   if (dir == 0) {
01584     if (j == m_nDiv[1] || k == m_nDiv[2]) return 0.; // return outside
01585     iEdge = k*m_nDiv[0]*m_nDiv[1] + j*m_nDiv[0] + i;
01586   }
01587   else if (dir == 1) {
01588     if (i == m_nDiv[0] || k == m_nDiv[2]) return 0.; // return outside
01589     iEdge = i*m_nDiv[1]*m_nDiv[2] + k*m_nDiv[1] + j;
01590   }
01591   else if (dir == 2) {
01592     if (i == m_nDiv[0] || j == m_nDiv[1]) return 0.; // return outside
01593     iEdge = j*m_nDiv[0]*m_nDiv[2] + i*m_nDiv[2] + k;
01594   }
01595   else return -1.;
01596 
01597   EdgeStatus edgeStatus = m_vnEdgeStatus[dir][iEdge];
01598   if (edgeStatus == INSIDE) return 1.;
01599   else if (edgeStatus == OUTSIDE) return 0.;
01600   else return -1.;
01601 }
01602 
01603 bool EBMesher::move_ray(int& nIntersect, double* startPnt, double* endPnt,
01604                       double tol, int dir, bool bSpecialCase)
01605 {
01606   //int i, nIteration;
01607   int i;
01608   bool bMove = true;
01609 
01610   if (bSpecialCase) m_iOverlap = 0;
01611 
01612   while (bMove) {
01613     if (m_nIteration > 10) {
01614       if (bSpecialCase) return true; // special case
01615       else if (m_bUseGeom) { // shared/overlapped, faceting case
01616         m_mhOverlappedSurf[m_iOverlap] = 0;
01617         m_mhOverlappedSurf[m_iOverlap + 1] = 0;
01618         return true;
01619       }
01620       else {
01621         return false;
01622       }
01623     }
01624 
01625     for (i = 0; i < 3; i++) {
01626       if (i != dir) {
01627         startPnt[i] += tol;
01628         endPnt[i] += tol;
01629       }
01630     }
01631 
01632     moab::CartVect ray(endPnt[0] - startPnt[0], endPnt[1] - startPnt[1], endPnt[2] - startPnt[2]);
01633     double rayLength = ray.length();
01634     moab::ErrorCode rVal;
01635     ray.normalize();
01636     m_vIntersection.clear();
01637     m_vhInterSurf.clear();
01638     m_vhInterFacet.clear();
01639     
01640     std::vector<double> temp_intersects;
01641     moab::OrientedBoxTreeTool::TrvStats stats;
01642     if (m_bUseGeom) {
01643       rVal = m_hObbTree->ray_intersect_sets(temp_intersects, m_vhInterSurf,
01644                                             m_vhInterFacet, m_hTreeRoot, tol,
01645                                             startPnt, ray.array(), &rayLength,
01646                                             &stats);
01647     }
01648     else { // facet input
01649       std::vector<moab::EntityHandle> dum_facets_out;
01650       rVal = m_hObbTree->ray_intersect_triangles(temp_intersects, dum_facets_out,
01651                                                  m_hTreeRoot, tol,
01652                                                  startPnt, ray.array(), &rayLength,
01653                                                  &stats);
01654       m_vhInterSurf.resize(temp_intersects.size());
01655     }
01656 
01657     if (moab::MB_SUCCESS != rVal) {
01658       std::cerr << "Failed : ray-triangle intersection." << std::endl;
01659       return false;
01660     }
01661 
01662     nIntersect = temp_intersects.size();
01663     m_vIntersection.resize(nIntersect);
01664     for (i = 0; i < nIntersect; i++) {
01665       IntersectDist temp_inter_dist(temp_intersects[i], i);
01666       m_vIntersection[i] = temp_inter_dist;
01667     }
01668     std::sort(m_vIntersection.begin(), m_vIntersection.end(), less_intersect);
01669     bMove = is_ray_move_and_set_overlap_surf(bSpecialCase);
01670     m_nIteration++;
01671   }
01672 
01673   std::cout << "ray is moved successfully." << std::endl;
01674 
01675   return true;
01676 }
01677 
01678 bool EBMesher::is_ray_move_and_set_overlap_surf(bool& bSpecialCase)
01679 {
01680   int nInter = m_vIntersection.size() - 1;
01681 
01682   while (m_iOverlap < nInter) {
01683     if (m_vIntersection[m_iOverlap + 1].distance - m_vIntersection[m_iOverlap].distance < 1e-7) {
01684       if (m_bUseGeom) {
01685         moab::EntityHandle h1 = m_vhInterSurf[m_vIntersection[m_iOverlap].index];
01686         moab::EntityHandle h2 = m_vhInterSurf[m_vIntersection[m_iOverlap + 1].index];
01687         
01688         if (h1 == h2) { // remove too close case
01689           bSpecialCase = false;
01690           return true;
01691         }
01692         else if (m_nIteration < 10) { // when ray intesect shared edge by 2 triangles
01693           bSpecialCase = true;
01694           return true;
01695         }
01696         else { // overlapped surfaces
01697           m_mhOverlappedSurf[h1] = 0;
01698           m_mhOverlappedSurf[h2] = 0;
01699           m_nIteration = 0;
01700           m_iOverlap++;
01701         }
01702       }
01703       else {
01704         if (m_nIteration < 10) { // when ray intesect shared edge by 2 triangles
01705           bSpecialCase = true;
01706           return true;
01707         }
01708         else { // overlapped surfaces
01709           //bSpecialCase = false;
01710           //m_nIteration = 0;
01711           //m_iOverlap++;
01712           m_vIntersection.erase(m_vIntersection.begin() + m_iOverlap + 1);
01713           nInter = m_vIntersection.size();
01714           m_vhInterSurf.resize(nInter);
01715           //m_nIteration = 0;
01716         }
01717       }
01718     }
01719     else m_iOverlap++;
01720   }
01721 
01722   return false;
01723 }
01724 
01725 void EBMesher::remove_intersection_duplicates()
01726 {
01727   int ind = 0;
01728   int nInter = m_vIntersection.size() - 1;
01729 
01730   while (ind < nInter) {
01731     if (m_vIntersection[ind + 1].distance - m_vIntersection[ind].distance < 1e-7) {
01732       moab::EntityHandle h1 = m_vhInterSurf[m_vIntersection[ind].index];
01733       moab::EntityHandle h2 = m_vhInterSurf[m_vIntersection[ind + 1].index];
01734       
01735       if (h1 != h2) { // overlapped/shared surfaces
01736         std::vector<iBase_EntitySetHandle> parents_out1, parents_out2;
01737         int err = mk_core()->imesh_instance()->getPrnts(reinterpret_cast<iBase_EntitySetHandle> (h1), 1, parents_out1);
01738         IBERRCHK(err, "Failed to get number of surface parents.\n");
01739         err = mk_core()->imesh_instance()->getPrnts(reinterpret_cast<iBase_EntitySetHandle> (h2), 1, parents_out2);
01740         IBERRCHK(err, "Failed to get number of surface parents.\n");
01741         if (parents_out1.size() == 1 && parents_out2.size() == 1) {
01742           if (parents_out1[0] != parents_out2[0]) { // parent volues are also different
01743             m_mhOverlappedSurf[h1] = 0;
01744             m_mhOverlappedSurf[h2] = 0;
01745           }
01746         }
01747       }
01748       m_vIntersection.erase(m_vIntersection.begin() + ind + 1);
01749       nInter--;
01750     }
01751     else ind++;
01752   }
01753 }
01754 
01755 bool EBMesher::get_volume_fraction(int vol_frac_div)
01756 {
01757   std::cout << "starting get_volume_fraction." << std::endl;
01758   int i, j, k, l, n, iHex, dir, nIntersect, rayIndex[3], index[3],
01759     otherDir1, otherDir2, err, nParent;
01760   (void) otherDir1;
01761   (void) otherDir2;
01762   double tolerance = 1e-12, dDistance, elem_origin[3],
01763     elem_interval_size[3], startPnt[3], endPnt[3], rayMaxEnd[3];
01764   moab::ErrorCode rval;
01765   bool bOverlapSecond, bClosed = true;
01766   std::vector<iBase_EntityHandle> edge_handles;
01767 
01768   double dTotEdgeElem = 0.;
01769   for (i = 0; i < 3; i++) {
01770     rayMaxEnd[i] = m_origin_coords[i] + m_nDiv[i]*m_dIntervalSize[i];
01771     dTotEdgeElem += m_dIntervalSize[i]*(vol_frac_div + 1)*(vol_frac_div + 1);
01772   }
01773 
01774   // get all hexes
01775   std::vector<iBase_EntityHandle> hex_handles;
01776   err = mk_core()->imesh_instance()->getEntities(m_hRootSet, iBase_REGION,
01777                                                  iMesh_HEXAHEDRON, hex_handles);
01778   IBERRCHK(err, "Failed to get hexes.\n");
01779 
01780   int hex_size = hex_handles.size();
01781   /*
01782   std::vector<int> hex_status(hex_size);
01783   //int* hex_status = NULL;
01784   err = mk_core()->imesh_instance()->getIntArrData(&hex_handles[0], hex_size,
01785                                                    m_elemStatusTag, &hex_status[0]);
01786   */
01787   int error;
01788   std::vector<int> frac_length;
01789   int* hex_status = new int[hex_size];
01790   int hex_status_alloc = sizeof(int)*hex_size;
01791   int hex_status_size = 0;
01792   iMesh_getIntArrData(m_mesh, &hex_handles[0],
01793                       hex_size, m_elemStatusTag,
01794                       &hex_status, &hex_status_alloc,
01795                       &hex_status_size, &error);
01796   ERRORRF("Failed to get hex status.\n");
01797 
01798   for (n = 0; n < hex_size; n++) {
01799     if (hex_status[n] == 2) { // just boundary
01800       std::map<moab::EntityHandle, VolFrac> vol_fraction;
01801       std::map<moab::EntityHandle, VolFrac>::iterator vf_iter;
01802       std::map<moab::EntityHandle, VolFrac>::iterator vf_end_iter;
01803       iHex = moab_instance()->id_from_handle(reinterpret_cast<moab::EntityHandle>
01804                                              (hex_handles[n])) - m_iStartHex;
01805       index[0] = (iHex%(m_nDiv[0]*m_nDiv[1]))%m_nDiv[0];
01806       index[1] = (iHex%(m_nDiv[0]*m_nDiv[1]))/m_nDiv[0];
01807       index[2] = iHex/m_nDiv[0]/m_nDiv[1];
01808       
01809       for (i = 0; i < 3; i++) {
01810         elem_origin[i] = m_origin_coords[i] + index[i]*m_dIntervalSize[i];
01811         elem_interval_size[i] = m_dIntervalSize[i]/vol_frac_div;
01812       }
01813       
01814       for (dir = 0; dir < 3; dir++) { // x, y, z directions
01815         otherDir1 = (dir + 1)%3;
01816         otherDir2 = (dir + 2)%3;
01817         
01818         for (j = 0; j < vol_frac_div + 1; j++) {
01819           for (i = 0; i < vol_frac_div + 1; i++) {
01820             // get ray start and end points
01821             if (dir == 0) { // x coordinate ray
01822               rayIndex[0] = 0;
01823               rayIndex[1] = i;
01824               rayIndex[2] = j;
01825             }
01826             else if (dir == 1) { // y coordinate ray
01827               rayIndex[0] = j;
01828               rayIndex[1] = 0;
01829               rayIndex[2] = i;
01830             }
01831             else if (dir == 2) { // z coordinate ray
01832               rayIndex[0] = i;
01833               rayIndex[1] = j;
01834               rayIndex[2] = 0;
01835             }
01836             else IBERRCHK(iBase_FAILURE, "Ray direction should be from 0 to 2.");
01837             
01838             for (k = 0; k < 3; k++) {
01839               if (k == dir) {
01840                 startPnt[k] = elem_origin[k];
01841                 endPnt[k] = startPnt[k] + m_dIntervalSize[k];
01842               }
01843               else {
01844                 startPnt[k] = elem_origin[k] + rayIndex[k]*elem_interval_size[k];
01845                 endPnt[k] = startPnt[k];
01846               }
01847             }
01848             
01849             // ray-tracing
01850             if (!fire_ray(nIntersect, startPnt, endPnt,
01851                           tolerance, dir, m_dIntervalSize[dir])) return iBase_FAILURE;
01852             
01853             if (nIntersect > 0) { // ray is intersected
01854               bOverlapSecond = false;
01855               bClosed = true;
01856               for (k = 0; k < nIntersect; k++) {
01857                 std::vector<moab::EntityHandle> parents;
01858                 //moab::Range parents;
01859                 rval = moab_instance()->get_parent_meshsets(m_vhInterSurf[m_vIntersection[k].index],
01860                                                             parents);
01861                 if (rval != moab::MB_SUCCESS) {
01862                   std::cerr << "Couldn't get parents." << std::endl;
01863                   return false;
01864                 }
01865 
01866                 nParent = parents.size();
01867                 dDistance = m_vIntersection[k].distance;
01868                 
01869                 if (is_shared_overlapped_surf(k)) {
01870                   if (nParent == 2) { // shared surface
01871                     for (l = 0; l < 2; l++) {
01872                       if (l == 1) {
01873                         dDistance *= -1.;
01874                         bClosed = false;
01875                       }
01876                       else bClosed = true;
01877 
01878                       vf_iter = vol_fraction.find(parents[l]);
01879                       if (vf_iter == vol_fraction.end()) {
01880                         VolFrac temp_vf(dDistance, bClosed);
01881                         vol_fraction[parents[l]] = temp_vf;
01882                         //std::cout << "iHex=" << iHex << ", vh="
01883                         //        << parents[l] << ", dir=" << dir << ", dDistance1="
01884                         //        << dDistance << ", vol="
01885                         //        << temp_vf.vol << std::endl;
01886                       }
01887                       else {
01888                         vf_iter->second.vol += dDistance;
01889                         vf_iter->second.closed = bClosed;
01890                         //std::cout << "iHex=" << iHex << ", vh="
01891                         //        << vf_iter->first << ", dir=" << dir << ", dDistance1="
01892                         //                         << dDistance << ", vol="
01893                         //        << vf_iter->second.vol << std::endl;
01894                       }
01895                     }
01896                   }
01897                   else if (nParent == 1) { // overlapped surface
01898                     if (bOverlapSecond) { // second intersection
01899                       /*
01900                       for (int m = 0; m < 3; m++) {
01901                         ePnt[m] = startPnt[m];
01902                       }
01903                       ePnt[dir] += dDistance;
01904                       */
01905                       dDistance *= -1.;
01906                       bClosed = false;
01907                       bOverlapSecond = false;
01908                     }
01909                     else {
01910                       /*
01911                       // make edge
01912                       for (int m = 0; m < 3; m++) {
01913                         if (bClosed) ePnt[m] = startPnt[m];
01914                         ePnt[m + 3] = ePnt[m];
01915                       }
01916                       ePnt[dir + 3] += dDistance;
01917                       if (!make_edge(ePnt, edge_handles)) return iBase_FAILURE;
01918                       */
01919                       bOverlapSecond = true;// first intersection
01920                       bClosed = true;
01921                     }
01922                     
01923                     vf_iter = vol_fraction.find(parents[0]);
01924                     if (vf_iter == vol_fraction.end()) {
01925                       VolFrac temp_vf(dDistance, bClosed);
01926                       vol_fraction[parents[0]] = temp_vf;
01927                       //std::cout << "iHex=" << iHex << ", vh="
01928                       //        << parents[0] << ", dDistance2="
01929                       //        << dDistance << ", vol="
01930                       //        << temp_vf.vol << std::endl;
01931                     }
01932                     else {
01933                       vf_iter->second.vol += dDistance;
01934                       vf_iter->second.closed = bClosed;
01935                       //std::cout << "iHex=" << iHex << ", vh="
01936                       //        << vf_iter->first << ", dir=" << dir << ", dDistance2="
01937                       //        << dDistance << ", vol="
01938                       //        << vf_iter->second.vol << std::endl;
01939                     }
01940                   }
01941                   else return iBase_FAILURE;
01942                 }
01943                 else { // normal case
01944                   if (nParent != 1) return iBase_FAILURE;
01945 
01946                   if (!is_same_direct_to_ray(k, dir)) {
01947                     dDistance *= -1.; // outside
01948                     bClosed = false;
01949                   }
01950                   else bClosed = true;
01951                   
01952                   vf_iter = vol_fraction.find(parents[0]);
01953                   if (vf_iter == vol_fraction.end()) {
01954                     VolFrac temp_vf(dDistance, bClosed);
01955                     vol_fraction[parents[0]] = temp_vf;
01956                     //std::cout << "iHex=" << iHex << ", vh="
01957                     //        << parents[0] << ", dir=" << dir << ", dDistance3="
01958                     //        << dDistance << ", vol="
01959                     //        << temp_vf.vol << std::endl;
01960                   }
01961                   else {
01962                     vf_iter->second.vol += dDistance;
01963                     vf_iter->second.closed = bClosed;
01964                     //std::cout << "iHex=" << iHex << ", vh="
01965                     //        << vf_iter->first << ", dir=" << dir << ", dDistance3="
01966                     //        << dDistance << ", vol="
01967                     //        << vf_iter->second.vol << std::endl;
01968                   }
01969                 }
01970               }
01971 
01972               // if fraction is not closed, add interval size to close
01973               vf_iter = vol_fraction.begin();
01974               vf_end_iter = vol_fraction.end();
01975               for (; vf_iter != vf_end_iter; vf_iter++) {
01976                 if (!vf_iter->second.closed) {
01977                   vf_iter->second.vol += m_dIntervalSize[dir];
01978                   vf_iter->second.closed = true;
01979                   /*
01980                   for (int m = 0; m < 3; m++) {
01981                     ePnt[m + 3] = startPnt[m];
01982                   }
01983                   ePnt[dir + 3] += m_dIntervalSize[dir];
01984                   if (!make_edge(ePnt, edge_handles)) return iBase_FAILURE;
01985                   */
01986                   //std::cout << "iHex=" << iHex << ", vh="
01987                   //        << vf_iter->first << ", dir=" << dir << ", dDistance4="
01988                   //        << m_dIntervalSize[dir] << ", vol="
01989                   //        << vf_iter->second.vol << std::endl;
01990                 }
01991               }
01992             }
01993             else { // decide if it is fully inside
01994               if (!fire_ray(nIntersect, startPnt, rayMaxEnd,
01995                             tolerance, dir, m_nDiv[dir]*m_dIntervalSize[dir])) return iBase_FAILURE;
01996               
01997               if (nIntersect > 0) { // fully inside
01998                 if (is_shared_overlapped_surf(0) || // shared or overlapped
01999                     is_same_direct_to_ray(0, dir)) { // other inside case
02000                   moab::Range parents;
02001                   rval = moab_instance()->get_parent_meshsets(m_vhInterSurf[m_vIntersection[0].index],
02002                                                               parents);
02003                   if (rval != moab::MB_SUCCESS) {
02004                     std::cerr << "Couldn't get parents." << std::endl;
02005                     return false;
02006                   }
02007 
02008                   moab::EntityHandle parent_entity = parents.pop_front();
02009                   vf_iter = vol_fraction.find(parent_entity);
02010                   if (vf_iter == vf_end_iter) {
02011                     VolFrac temp_vf(m_dIntervalSize[dir], true);
02012                     vol_fraction[parent_entity] = temp_vf;
02013                     //std::cout << "iHex=" << iHex << ", vh="
02014                     //        << parents[0] << ", dir=" << dir << ", dDistance5="
02015                     //                             << dDistance << ", vol="
02016                     //        << temp_vf.vol << std::endl;
02017                   }
02018                   else {
02019                     vf_iter->second.vol += m_dIntervalSize[dir];
02020                     vf_iter->second.closed = bClosed;
02021                     //std::cout << "iHex=" << iHex << ", vh="
02022                     //        << vf_iter->first << ", dir=" << dir << ", dDistance5="
02023                     //                         << dDistance << ", vol="
02024                     //        << vf_iter->second.vol << std::endl;
02025                   }
02026                 }
02027               }
02028             }
02029           }
02030         }
02031       }
02032 
02033       int vol_frac_length = vol_fraction.size();
02034       std::vector<int> vol_frac_id(vol_frac_length);
02035       std::vector<double> vol_frac(vol_frac_length);
02036       int vol_frac_id_size = vol_frac_length*sizeof(int);
02037       int vol_frac_size = vol_frac_length*sizeof(double);
02038       vf_iter = vol_fraction.begin();
02039       vf_end_iter = vol_fraction.end();
02040       for (int m = 0; vf_iter != vf_end_iter; vf_iter++, m++) {
02041         vol_frac_id[m] = vf_iter->first;
02042         vol_frac[m] = vf_iter->second.vol/dTotEdgeElem;
02043         //std::cout << "iHex=" << iHex << ", i=" << index[0]
02044         //        << ", j=" << index[1] << ", k=" << index[2]
02045         //        << ", vol_handle=" << vf_iter->first
02046         //        << ", vol=" << vf_iter->second.vol
02047         //        << ", vol_fraction=" << vf_iter->second.vol/dTotEdgeElem
02048         //        << std::endl;
02049       }
02050       const void* vol_frac_ids = &vol_frac_id[0];
02051       const void* vol_fracs = &vol_frac[0];
02052 
02053       // set volume fraction information as tag
02054       rval = moab_instance()->tag_set_data(reinterpret_cast<moab::Tag> (m_volFracLengthTag),
02055                                            reinterpret_cast<moab::EntityHandle*> (&hex_handles[n]),
02056                                            1, &vol_frac_length);
02057       MBERRCHK(rval, mk_core()->moab_instance());
02058 
02059       rval = moab_instance()->tag_set_by_ptr(reinterpret_cast<moab::Tag> (m_volFracHandleTag),
02060                                            reinterpret_cast<moab::EntityHandle*> (&hex_handles[n]),
02061                                            1, &vol_frac_ids, &vol_frac_id_size);
02062       MBERRCHK(rval, mk_core()->moab_instance());
02063 
02064       rval = moab_instance()->tag_set_by_ptr(reinterpret_cast<moab::Tag> (m_volFracTag),
02065                                            reinterpret_cast<moab::EntityHandle*> (&hex_handles[n]),
02066                                            1, &vol_fracs, &vol_frac_size);
02067       MBERRCHK(rval, mk_core()->moab_instance());
02068 
02069       if (debug_ebmesh) {
02070         for (int v = 0; v < vol_frac_length; v++) {
02071           std::cout << "#_boundary_hex_handles=" << hex_handles[n]
02072                     << ",vol_frac_handle[" << v << "]=" << vol_frac_id[v]
02073                     << ",vol_fracs[" << v << "]=" << vol_frac[v]
02074                     << std::endl;
02075         }
02076       }
02077     }
02078   }  
02079 
02080   std::cout << "end get_volume_fraction." << std::endl;
02081   delete [] hex_status;
02082 
02083   return true;
02084 }
02085 
02086 bool EBMesher::is_same_direct_to_ray(int i, int dir)
02087 {
02088   moab::CartVect coords[3], normal(0.0), ray_dir(rayDir[dir]);
02089   const moab::EntityHandle *conn;
02090   int len;
02091 
02092   // get triangle normal vector
02093   moab::ErrorCode rval = moab_instance()->get_connectivity(m_vhInterFacet[m_vIntersection[i].index],
02094                                                            conn, len);
02095   assert(moab::MB_SUCCESS == rval && 3 == len);
02096   
02097   rval = moab_instance()->get_coords(conn, 3, coords[0].array());
02098   (void) rval;
02099   assert(moab::MB_SUCCESS == rval);
02100 
02101   coords[1] -= coords[0];
02102   coords[2] -= coords[0];
02103   normal += coords[1] * coords[2];
02104   normal.normalize();
02105 
02106   return angle(normal, ray_dir) < .5*PI;
02107 }
02108 
02109 void EBMesher::set_obb_tree_box_dimension()
02110 {
02111   std::cout << "starting to construct obb tree." << std::endl;
02112   moab::ErrorCode rval = m_GeomTopoTool->construct_obb_trees(m_bUseWholeGeom);
02113   MBERRCHK(rval, mk_core()->moab_instance());
02114 
02115   m_hObbTree = m_GeomTopoTool->obb_tree();
02116   m_hTreeRoot = m_GeomTopoTool->get_one_vol_root();
02117 
02118   moab::CartVect box_center, box_axis1, box_axis2, box_axis3;
02119   for (int i = 0; i < 3; i++) {
02120     m_minCoord[i] = HUGE_VAL;
02121     m_maxCoord[i] = 0.;
02122   }
02123   rval = m_hObbTree->box(m_hTreeRoot, box_center.array(),
02124                          box_axis1.array(), box_axis2.array(),
02125                          box_axis3.array());
02126   MBERRCHK(rval, mk_core()->moab_instance());
02127   
02128   // cartesian box corners
02129   moab::CartVect corners[8] = {box_center + box_axis1 + box_axis2 + box_axis3,
02130                                box_center + box_axis1 + box_axis2 - box_axis3,
02131                                box_center + box_axis1 - box_axis2 + box_axis3,
02132                                box_center + box_axis1 - box_axis2 - box_axis3,
02133                                box_center - box_axis1 + box_axis2 + box_axis3,
02134                                box_center - box_axis1 + box_axis2 - box_axis3,
02135                                box_center - box_axis1 - box_axis2 + box_axis3,
02136                                box_center - box_axis1 - box_axis2 - box_axis3};
02137   
02138   // get the max, min cartesian box corners
02139   for (int i = 0; i < 8; i++) {
02140     for (int j = 0; j < 3; j++) {
02141       if (corners[i][j] < m_minCoord[j]) m_minCoord[j] = corners[i][j];
02142       if (corners[i][j] > m_maxCoord[j]) m_maxCoord[j] = corners[i][j];
02143     }
02144   }
02145   std::cout << "finished to construct obb tree." << std::endl;
02146 }
02147 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines