MeshKit  1.0
CAMALTetMesher.cpp
Go to the documentation of this file.
00001 #include "meshkit/CAMALTetMesher.hpp"
00002 #include "meshkit/MeshScheme.hpp"
00003 #include "meshkit/ModelEnt.hpp"
00004 #include "moab/Interface.hpp"
00005 #include "moab/ReadUtilIface.hpp"
00006 #include "moab/Range.hpp"
00007 #include "meshkit/RegisterMeshOp.hpp"
00008 #include "CAMALSurfEval.hpp"
00009 #include "CAMALSizeEval.hpp"
00010 #include "CMLTetMesher.hpp"
00011 #include "RefEntity.hpp"
00012 #include <vector>
00013 
00014 #ifdef PARALLEL
00015 #ifdef HAVE_PARALLEL_MOAB
00016 #include "moab/ParallelComm.hpp"
00017 #endif
00018 #endif
00019 
00020 const bool debug_camaltet = false;
00021 
00022 namespace MeshKit
00023 {
00024 
00025 moab::EntityType CAMALTetMesher::meshTps[] = {moab::MBVERTEX, moab::MBTET, moab::MBMAXTYPE};
00026 
00027 CAMALTetMesher::CAMALTetMesher(MKCore *mk_core, const MEntVector &me_vec)
00028         : MeshScheme(mk_core, me_vec)
00029 {
00030 }
00031 
00032 CAMALTetMesher::~CAMALTetMesher()
00033 {
00034 }
00035 
00036 MeshOp *CAMALTetMesher::get_tri_mesher() 
00037 {
00038   MeshOpProxy * mproxy = mk_core()->meshop_proxy("CAMALTriAdvance") ;
00039   return mk_core()->construct_meshop(mproxy);
00040 }
00041 
00042 void CAMALTetMesher::setup_this()
00043 {
00044   MeshOp *tri_mesher = NULL;
00045   std::vector<MeshOp*> meshops;
00046   MEntVector surfs;
00047   
00048   for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
00049     ModelEnt *me = (*sit).first;
00050     if (me->dimension() != 3) throw Error(MK_BAD_INPUT, "Tet mesher assigned to an entity with dimension != 3.");
00051     
00052       // get the bounding faces
00053     surfs.clear();
00054     me->get_adjacencies(2, surfs);
00055     
00056       // check the mesh scheme on each one; if there is one, verify it can generate tris; if there
00057       // isn't one, make one
00058     bool inserted = false;
00059 
00060     for (MEntVector::iterator vit = surfs.begin(); vit != surfs.end(); vit++) {
00061       if ((*vit)->is_meshops_list_empty()) {
00062           // get a tri mesher if we haven't already
00063         if (!tri_mesher) tri_mesher = get_tri_mesher();
00064         
00065           // add this surface to it, and if first for the volume, make sure it's added upstream
00066         tri_mesher->add_modelent(*vit);
00067         if (!inserted) {
00068           mk_core()->insert_node(tri_mesher, this);
00069           inserted = true;
00070         }
00071       } // if no meshops
00072     } // over surfs
00073   } // over vols
00074 }
00075 
00076 void CAMALTetMesher::execute_this()
00077 {
00078   for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
00079       // make a me, for convenience
00080     ModelEnt *me = (*sit).first;
00081 
00082     // assemble bounding mesh
00083     std::vector<moab::EntityHandle> bdy;
00084     std::vector<int> group_sizes, senses, bdy_ids;
00085     me->boundary(2, bdy, &senses, &group_sizes);
00086     
00087       // get connectivity 
00088 
00089       // convert the handles to integers for input to CMLTetMesher
00090     moab::Range bdy_vrange;
00091     std::vector<double> coords;
00092     me->get_indexed_connect_coords(bdy, &senses, NULL, bdy_ids, coords, &bdy_vrange);
00093 
00094     CMLTetMesher tet_mesher;
00095     bool success = tet_mesher.set_boundary_mesh(bdy_vrange.size(), &coords[0], bdy.size(), &bdy_ids[0]);
00096     if (!success)
00097       ECERRCHK(MK_FAILURE, "Failed setting boundary mesh");
00098 
00099       // generate the mesh
00100     int num_tets, num_pts;
00101     success = tet_mesher.generate_mesh(num_pts, num_tets);
00102     if (!success) {
00103       if (debug_camaltet) print_debug(me, coords, bdy_vrange, bdy, group_sizes, bdy_ids);
00104       ECERRCHK(MK_FAILURE, "Trouble generating tet mesh.");
00105     }
00106 
00107     moab::Range &new_ents = (*sit).second;
00108     moab::ErrorCode rval;
00109 
00110       // resize the coords array, then get the coords of the new points
00111     assert(num_pts >= (int)bdy_vrange.size());
00112     if (num_pts > (int)bdy_vrange.size()) {
00113       coords.resize(3*(num_pts-bdy_vrange.size()));
00114       int pts_returned = tet_mesher.get_points_buf(coords.size(), &coords[0], bdy_vrange.size());
00115       if (pts_returned != num_pts-(int)bdy_vrange.size())
00116         ECERRCHK(MK_FAILURE, "Number of new points returned from TetMesher doesn't agree with previous value output.");
00117     
00118         // create the new vertices' entities 
00119       rval = mk_core()->moab_instance()->create_vertices(&coords[0], pts_returned, new_ents);
00120       MBERRCHK(rval, mk_core()->moab_instance());
00121     }
00122 
00123       // for tets, pre-allocate connectivity
00124     moab::ReadUtilIface *iface;
00125     rval = mk_core()-> moab_instance() -> query_interface(iface);
00126     MBERRCHK(rval, mk_core()->moab_instance());         
00127 
00128       //create the tris, get a direct ptr to connectivity
00129     moab::EntityHandle starth, *connect;// *tmp_connect;
00130     rval = iface->get_element_connect(num_tets, 4, moab::MBTET, 1, starth, connect);
00131     MBERRCHK(rval, mk_core()->moab_instance());
00132 
00133       // read connectivity directly into that array, as int's
00134     int *connecti = (int*) connect;
00135     int tets_returned = tet_mesher.get_tets_buf(4*num_tets, connecti);
00136     if (tets_returned != num_tets)
00137       ECERRCHK(MK_FAILURE, "Number of new tets returned from TetMesher doesn't agree with previous value output.");
00138 
00139       // put vertex handles into an indexible array
00140     std::vector<moab::EntityHandle> bdy_vvec;
00141     std::copy(bdy_vrange.begin(), bdy_vrange.end(), std::back_inserter(bdy_vvec));
00142     std::copy(new_ents.begin(), new_ents.end(), std::back_inserter(bdy_vvec));
00143     
00144       // now convert vertex indices into handles in-place, working from the back
00145     for (int i = 4*num_tets-1; i >= 0; i--) {
00146       assert(connecti[i] >= 0 && connecti[i] < (int)bdy_vvec.size());
00147       connect[i] = bdy_vvec[connecti[i]];
00148     }
00149     
00150       // put new tris into new entity range, then commit the mesh
00151     new_ents.insert(starth, starth+num_tets-1);
00152     me->commit_mesh(new_ents, COMPLETE_MESH);
00153   }
00154 }
00155 
00156 void CAMALTetMesher::print_debug(ModelEnt *me, std::vector<double> &coords,
00157                                  moab::Range &bdy_vrange, std::vector<moab::EntityHandle> &bdy,
00158                                  std::vector<int> &group_sizes, std::vector<int> &bdy_ids)
00159 {
00160   std::cout << "Volume_boundary_mesh: mesh_size = "
00161             << me->mesh_interval_size() << std::endl;
00162   
00163   for (int i = 0; i < (int)bdy_vrange.size(); i++) {
00164     std::cout << coords[3 * i] << "  " << coords[3 * i + 1]
00165               << "  " << coords[3 * i + 2] << std::endl;
00166   }
00167   
00168   std::cout << "bdy_vertex_size:" << bdy_vrange.size()
00169             << ", group_size:" << group_sizes.size() << std::endl;
00170   
00171   int index = 0;
00172   for (size_t i = 0; i < group_sizes.size(); i++) {
00173     int g_size = group_sizes[i];
00174     std::cout << "boundary_order_group" << i + 1 << ", group_size="
00175               << g_size << std::endl;
00176     for (int j = 0; j < g_size; j++) {
00177       std::cout << bdy_ids[index + j] << " ";
00178     }
00179     std::cout << std::endl;
00180     index += g_size;
00181   }
00182   
00183   moab::ErrorCode rval;
00184   moab::EntityHandle outset;
00185   std::string outfile;
00186   std::stringstream ss;
00187   RefEntity* entity = reinterpret_cast<RefEntity*> (me->geom_handle());
00188   ss << "CAMALTet_boundary";
00189   ss << entity->id();
00190   ss << "_";
00191 #ifdef PARALLEL
00192 #ifdef HAVE_PARALLEL_MOAB
00193   moab::ParallelComm* pcomm = moab::ParallelComm::get_pcomm(mk_core()->moab_instance(), 0);
00194   ss << "proc";
00195   ss << pcomm->proc_config().proc_rank();
00196 #endif
00197 #endif
00198   ss >> outfile;
00199   outfile += ".vtk";
00200   rval = mk_core()->moab_instance()->create_meshset(0, outset);
00201   MBERRCHK(rval, mk_core()->moab_instance());
00202   rval = mk_core()->moab_instance()->add_entities(outset, &bdy[0], bdy.size());
00203   MBERRCHK(rval, mk_core()->moab_instance());
00204   rval = mk_core()->moab_instance()->write_mesh(outfile.c_str(), &outset, 1);
00205   MBERRCHK(rval, mk_core()->moab_instance());
00206   std::cout << outfile.c_str() << " is saved." << std::endl;
00207 }
00208 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines