MeshKit  1.0
SCDMesh.cpp
Go to the documentation of this file.
00001 //-----------------------------------C++-------------------------------------//
00002 // File: algs/SCDMesh.cpp
00003 // Author: Stuart R. Slattery
00004 // Wednesday February 2 16:15:16 2011
00005 // Brief: SCDMesh member definitions
00006 //---------------------------------------------------------------------------//
00007 
00008 #include "meshkit/SCDMesh.hpp"
00009 #include "meshkit/MKCore.hpp"
00010 #include "meshkit/ModelEnt.hpp"
00011 #include "meshkit/RegisterMeshOp.hpp"
00012 
00013 #include <vector>
00014 
00015 #include "moab/Core.hpp"
00016 #include "moab/Range.hpp"
00017 #include "moab/Interface.hpp"
00018 #include "moab/ScdInterface.hpp"
00019 #include "moab/OrientedBoxTreeTool.hpp"
00020 #include "moab/CartVect.hpp"
00021 
00022 namespace MeshKit
00023 {
00024 //---------------------------------------------------------------------------//
00025 // Initialize the entity types for SCDMesh
00026   moab::EntityType SCDMesh_tps[] = {moab::MBVERTEX, moab::MBHEX, moab::MBMAXTYPE};
00027   const moab::EntityType* SCDMesh::output_types() 
00028     { return SCDMesh_tps; }
00029 
00030 //---------------------------------------------------------------------------//
00031 // Constructor for SCDMesh
00032   SCDMesh::SCDMesh(MKCore *mk_core, const MEntVector &me_vec)
00033     : MeshScheme(mk_core, me_vec)
00034   {
00035     boxIncrease = 0.0;
00036     useMeshGeom = false;
00037 
00038     // set bounding box size tag
00039     double bb_default[6] = { 0., 0., 0., 0., 0., 0. };
00040     rval = mk_core->moab_instance()->tag_get_handle("BOUNDING_BOX_SIZE",
00041                                                 6, moab::MB_TYPE_DOUBLE,
00042             bb_tag, moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT, bb_default);
00043     if (moab::MB_SUCCESS != rval && moab::MB_ALREADY_ALLOCATED != rval) 
00044       MBERRCHK(rval, mk_core->moab_instance());
00045   }
00046 
00047 //---------------------------------------------------------------------------//
00048 // setup function
00049 void SCDMesh::setup_this()
00050 {
00051   // if cfMesh chosen, check the grid entries to make sure they're correct
00052   if (gridType == 0) {
00053     if (coarse_i != fine_i.size()) {
00054       throw Error(MK_INCOMPLETE_MESH_SPECIFICATION, 
00055       "Number of coarse i divisions not equal to the number of fine i divisions.");
00056     }
00057     if (coarse_j != fine_j.size()) {
00058       throw Error(MK_INCOMPLETE_MESH_SPECIFICATION, 
00059       "Number of coarse j divisions not equal to the number of fine j divisions.");
00060       }
00061     if (coarse_k != fine_k.size()) {
00062       throw Error(MK_INCOMPLETE_MESH_SPECIFICATION, 
00063       "Number of coarse k divisions not equal to the number of fine k divisions.");
00064     }
00065   }
00066 
00067   // do setup for the model entities
00068   for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++) {
00069     ModelEnt *me = mit->first;
00070 
00071     // check if the entity is already meshed
00072     if (me->get_meshed_state() >= COMPLETE_MESH || me->mesh_intervals() > 0) continue;
00073   }
00074 }
00075 
00076 //---------------------------------------------------------------------------//
00077 // Structured cartesian mesh generation execution
00078   void SCDMesh::execute_this()
00079   {
00080     for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++) {
00081       ModelEnt *me = mit->first;
00082 
00083       /* Generate the bounding box */
00084 
00085       // Cartesian case
00086       if (axisType == 0) {
00087         if (geomType == 0) {
00088           set_cart_box_all();
00089         }
00090         else if (geomType == 1) {
00091           set_cart_box_individual(me);
00092         }
00093         create_cart_edges();
00094       }
00095 
00096       /* Generate the mesh */
00097 
00098       // create mesh vertex coordinates
00099       create_vertex_coords();
00100 
00101       // full representation case
00102       if (interfaceType == 0) {
00103         create_full_mesh(mit->second);
00104       }
00105 
00106       // lighweight ScdInterface case
00107       if (interfaceType == 1) {
00108         create_light_mesh(mit->second);
00109       }
00110 
00111       // stop if mesh is created for whole geometry
00112       if (geomType == 0) break;
00113 
00114       // commit the mesh when finished
00115       me->commit_mesh(mit->second, COMPLETE_MESH);
00116     }
00117   }
00118 
00119 //---------------------------------------------------------------------------//
00120 // set the Cartesian bounding box dimension for the entire geometry
00121   void SCDMesh::set_cart_box_all()
00122   {
00123     if (!useMeshGeom) {
00124       gerr = mk_core()->igeom_instance()->getBoundBox(minCoord[0], minCoord[1], minCoord[2],
00125                                                       maxCoord[0], maxCoord[1], maxCoord[2]);
00126       IBERRCHK(gerr, "ERROR: couldn't get geometry bounding box");
00127     }
00128 
00129     // increase box size
00130     double box_size;
00131     for (int i = 0; i < 3; i++) {
00132       box_size = maxCoord[i] - minCoord[i];
00133       minCoord[i] -= box_size*boxIncrease; 
00134       maxCoord[i] += box_size*boxIncrease; 
00135     }
00136   }
00137   
00138 //---------------------------------------------------------------------------//
00139 // set the Cartesian bounding box dimension for an individual volume
00140   void SCDMesh::set_cart_box_individual(ModelEnt *this_me)
00141   {
00142     gerr = this_me->igeom_instance()->getBoundBox(minCoord[0], minCoord[1], minCoord[2],
00143                                                   maxCoord[0], maxCoord[1], maxCoord[2]);
00144     IBERRCHK(gerr, "ERROR: couldn't get geometry bounding box");
00145 
00146     // increase box size
00147     double box_size;
00148     for (int i = 0; i < 3; i++) {
00149       box_size = maxCoord[i] - minCoord[i];
00150       minCoord[i] -= box_size*boxIncrease; 
00151       maxCoord[i] += box_size*boxIncrease; 
00152     }
00153 
00154     // set bounding box size as tag
00155     moab::EntityHandle meshset = this_me->mesh_handle();
00156     double box_min_max[6] = { minCoord[0], minCoord[1], minCoord[2],
00157                               maxCoord[0], maxCoord[1], maxCoord[2] };
00158     rval = mk_core()->moab_instance()->tag_set_data(bb_tag, &meshset, 1,
00159                                                     box_min_max);
00160   }
00161 
00162 //---------------------------------------------------------------------------//
00163 // set the Cartesian bounding box dimension
00164 void SCDMesh::set_cart_box_min_max(double* min, double* max, double box_increase)
00165 {
00166   for (int i = 0; i < 3; i++) {
00167     double box_size = max[i] - min[i];
00168     minCoord[i] = min[i] - box_size*box_increase;
00169     maxCoord[i] = max[i] + box_size*box_increase;
00170   }
00171 }
00172 
00173 //---------------------------------------------------------------------------//
00174 // get the axis cartesian bounding box edges if cfmesh has been chosen
00175 void SCDMesh::create_cart_edges()
00176 {
00177   // generate the vector arrays for cartesian grid.
00178   i_arr.push_back(minCoord[0]);
00179   double icrs_size = (maxCoord[0] - minCoord[0]) / coarse_i;
00180   double ifn_size;
00181   unsigned int crsi;
00182   for (crsi = 0; crsi != coarse_i; crsi++) {
00183     ifn_size = icrs_size / fine_i[crsi];
00184     for (int fni = 1; fni != fine_i[crsi] + 1; fni++) {
00185       i_arr.push_back(minCoord[0] + crsi*icrs_size + fni*ifn_size);
00186     }
00187   }
00188   
00189   j_arr.push_back(minCoord[1]);
00190   double jcrs_size = (maxCoord[1] - minCoord[1]) / coarse_j;
00191   double jfn_size;
00192   unsigned int crsj;
00193   for (crsj = 0; crsj != coarse_j; crsj++) {
00194     jfn_size = jcrs_size / fine_j[crsj];
00195     for (int fnj = 1; fnj != fine_j[crsj] + 1; fnj++) {
00196       j_arr.push_back(minCoord[1] + crsj*jcrs_size + fnj*jfn_size);
00197     }
00198   }
00199   
00200   k_arr.push_back(minCoord[2]);
00201   double kcrs_size = (maxCoord[2] - minCoord[2]) / coarse_k;
00202   double kfn_size;
00203   unsigned int crsk;
00204   for (crsk = 0; crsk != coarse_k; crsk++) {
00205     kfn_size = kcrs_size / fine_k[crsk];
00206     for (int fnk = 1; fnk != fine_k[crsk] + 1; fnk++) {
00207       k_arr.push_back(minCoord[2] + crsk*kcrs_size + fnk*kfn_size);
00208     }
00209   }
00210 }
00211 
00212 //---------------------------------------------------------------------------//
00213 // create the mesh vertex coordinates
00214   void SCDMesh::create_vertex_coords()
00215   {
00216     num_i = i_arr.size();
00217     num_j = j_arr.size();
00218     num_k = k_arr.size();
00219     num_verts = num_i*num_j*num_k;
00220 
00221     full_coords.resize(3*num_verts);
00222     unsigned int idx;
00223 
00224     // generate the vertex coordinate array in an interleaved fashion
00225     for (int kval = 0; kval != num_k; kval++) {
00226       for (int jval = 0; jval != num_j; jval++) {
00227         for (int ival = 0; ival != num_i; ival++) {
00228           idx = ival + num_i*jval + num_i*num_j*kval;
00229           full_coords[3*idx] = i_arr[ival];
00230           full_coords[3*idx + 1] = j_arr[jval];
00231           full_coords[3*idx + 2] = k_arr[kval];
00232         }
00233       }
00234     }
00235   }
00236 
00237 //---------------------------------------------------------------------------//
00238 // create the full mesh representation
00239   void SCDMesh::create_full_mesh(moab::Range& me_range)
00240   {
00241     // create the vertices from the coordinates array
00242     moab::Range vtx_range;
00243     rval = mk_core()->moab_instance()->create_vertices(&full_coords[0],
00244                                                        num_verts,
00245                                                        vtx_range);
00246     MBERRCHK(rval, mk_core()->moab_instance());
00247 
00248     if (geomType == 0) me_range.merge(vtx_range);
00249 
00250     // generate hex entities from the vertices
00251     // the entity handles here should be contiguous with the i blocks 
00252     // moving the fastest and the k blocks moving the slowest
00253     moab::EntityHandle local_conn[8];
00254     unsigned int idx;
00255     for (int kv = 0; kv != num_k - 1; kv++) {
00256       for (int jv = 0; jv != num_j - 1; jv++) {
00257         for (int iv = 0; iv != num_i - 1; iv++) {
00258           idx = iv + jv*num_i + kv*num_i*num_j;
00259           local_conn[0] = vtx_range[ idx];
00260           local_conn[1] = vtx_range[ idx + 1];
00261           local_conn[2] = vtx_range[ idx + 1 + num_i];
00262           local_conn[3] = vtx_range[ idx +     num_i];
00263           local_conn[4] = vtx_range[ idx +             num_i*num_j];
00264           local_conn[5] = vtx_range[ idx + 1 +         num_i*num_j];
00265           local_conn[6] = vtx_range[ idx + 1 + num_i + num_i*num_j];
00266           local_conn[7] = vtx_range[ idx +     num_i + num_i*num_j];
00267 
00268           moab::EntityHandle this_elem;
00269           rval = mk_core()->moab_instance()->create_element( moab::MBHEX,
00270                                                              local_conn,
00271                                                              8,
00272                                                              this_elem);
00273           MBERRCHK(rval, mk_core()->moab_instance());
00274           if (geomType == 0) me_range.insert(this_elem);
00275         }
00276       }
00277     }
00278   }
00279 
00280 //---------------------------------------------------------------------------//
00281 // create mesh using light-weight MOAB ScdInterface
00282   void SCDMesh::create_light_mesh(moab::Range& me_range)
00283   {
00284     // create an instance of the ScdInterface
00285     rval = mk_core()->moab_instance()->query_interface(scdIface);
00286     MBERRCHK(rval, mk_core()->moab_instance());
00287 
00288     // create an instance of the ScdBox class
00289     // for now the coordinate indexing starts at 0
00290     // this should be changed to account for multiple geometric instances being meshed
00291     moab::ScdBox *scd_box;
00292     rval = scdIface->construct_box(moab::HomCoord(0, 0, 0, 1),
00293                                    moab::HomCoord(num_i-1, num_j-1, num_k-1, 1),
00294                                    &full_coords[0], 
00295                                    num_verts,
00296                                    scd_box);
00297     MBERRCHK(rval, mk_core()->moab_instance());
00298 
00299     // get rid of ScdInterface once we are done
00300     rval = mk_core()->moab_instance()->release_interface(scdIface);
00301     MBERRCHK(rval, mk_core()->moab_instance());
00302 
00303     // add created vertex and hexes
00304     if (geomType == 1) {
00305       moab::EntityHandle start;
00306       start = scd_box->start_vertex();
00307       moab::Range vert_range(start, start + scd_box->num_vertices()-1);
00308       me_range.merge(vert_range);
00309       start = scd_box->start_element();
00310       moab::Range hex_range(start, start + scd_box->num_elements()-1);
00311       me_range.merge(hex_range);
00312     }
00313   }
00314 
00315 //---------------------------------------------------------------------------//
00316 
00317 } // end namespace MeshKit
00318 
00319 //---------------------------------------------------------------------------//
00320 // end algs/SCDMesh.cpp 
00321 //---------------------------------------------------------------------------//
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines