MeshKit  1.0
CAMALTriAdvance.cpp
Go to the documentation of this file.
00001 #include "meshkit/CAMALTriAdvance.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 "CMLTriAdvance.hpp"
00011 #include "RefEntity.hpp"
00012 
00013 #ifdef PARALLEL
00014 #ifdef HAVE_PARALLEL_MOAB
00015 #include "moab/ParallelComm.hpp"
00016 #endif
00017 #endif
00018 
00019 #include <vector>
00020 
00021 const bool debug_camaltriadv = false;
00022 
00023 namespace MeshKit
00024 {
00025 
00026 moab::EntityType CAMALTriAdvance::meshTps[] = {moab::MBVERTEX, moab::MBTRI, moab::MBMAXTYPE};
00027 
00028 CAMALTriAdvance::CAMALTriAdvance(MKCore *mk_core, const MEntVector &me_vec)
00029         : MeshScheme(mk_core, me_vec)
00030 {
00031 }
00032 
00033 CAMALTriAdvance::~CAMALTriAdvance()
00034 {
00035 }
00036 
00037 void CAMALTriAdvance::setup_this()
00038 {
00039     // just call setup_boundary, since that's the only constraint we have
00040   setup_boundary();
00041 }
00042 
00043 void CAMALTriAdvance::execute_this()
00044 {
00045   for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
00046       // make a me, for convenience
00047     ModelEnt *me = (*sit).first;
00048 
00049       // create a surface evaluator for this modelent, and a size evaluator
00050     CAMALSurfEval cse(me);
00051 
00052     SizingFunction * sz = mk_core()->sizing_function(me->sizing_function_index());
00053     if (!sz->variable())
00054       sz = NULL; // no function variable
00055 
00056     CAMALSizeEval mesize(me->mesh_interval_size(), sz);
00057       // make sure the size isn't negative
00058     if (mesize.get_size() == -1) mesize.set_size(1.0);
00059     
00060     // assemble bounding mesh
00061     std::vector<moab::EntityHandle> bdy;
00062     std::vector<int> group_sizes, bdy_ids;
00063     me->boundary(0, bdy, NULL, &group_sizes);
00064 
00065       // convert the handles to integers for input to TriAdv
00066     moab::Range bdy_vrange;
00067     std::vector<double> coords;
00068     me->get_indexed_connect_coords(bdy, NULL, NULL, bdy_ids, coords, &bdy_vrange);
00069 
00070     // now construct the CAMAL mesher, and pass it initial conditions
00071     CMLTriAdvance triadv(&cse, &mesize);
00072     bool success = triadv.set_boundary_mesh(bdy_vrange.size(), &coords[0], group_sizes.size(), &group_sizes[0], &bdy_ids[0]);
00073     if (!success)
00074       ECERRCHK(MK_FAILURE, "Trouble setting boundary mesh.");
00075 
00076       // ok, now generate the mesh
00077     int num_pts, num_tris;
00078     success = triadv.generate_mesh(num_pts, num_tris);
00079     if (!success) {
00080       if (debug_camaltriadv) print_debug(me, coords, bdy_vrange, group_sizes, bdy_ids);
00081       ECERRCHK(MK_FAILURE, "Trouble generating tri mesh.");
00082     }
00083 
00084     moab::Range &new_ents = (*sit).second;
00085     moab::ErrorCode rval;
00086 
00087       // resize the coords array, then get the coords of the new points
00088     assert(num_pts >= (int)bdy_vrange.size());
00089     if (num_pts > (int)bdy_vrange.size()) {
00090       coords.resize(3*(num_pts-bdy_vrange.size()));
00091       unsigned int pts_returned = triadv.get_points_buf(coords.size(), &coords[0], bdy_vrange.size());
00092       if (pts_returned != num_pts-bdy_vrange.size()) 
00093         ECERRCHK(MK_FAILURE, "Number of new points returned from TriAdv doesn't agree with previous value output.");
00094     
00095         // create the new vertices' entities on the face
00096       rval = mk_core()->moab_instance()->create_vertices(&coords[0], pts_returned, new_ents);
00097       MBERRCHK(rval, mk_core()->moab_instance());
00098     }
00099 
00100       // for tris, pre-allocate connectivity
00101     moab::ReadUtilIface *iface;
00102     rval = mk_core()-> moab_instance() -> query_interface(iface);
00103     MBERRCHK(rval, mk_core()->moab_instance());         
00104 
00105       //create the tris, get a direct ptr to connectivity
00106     moab::EntityHandle starth, *connect;
00107     rval = iface->get_element_connect(num_tris, 3, moab::MBTRI, 1, starth, connect);
00108     MBERRCHK(rval, mk_core()->moab_instance());
00109 
00110       // read connectivity directly into that array, as int's
00111     int *connecti = (int*) connect;
00112     int tris_returned = triadv.get_tris_buf(3*num_tris, connecti);
00113     if (tris_returned != num_tris)
00114       ECERRCHK(MK_FAILURE, "Number of new tris returned from TriAdv doesn't agree with previous value output.");
00115 
00116       // put vertex handles into an indexible array
00117     std::vector<moab::EntityHandle> bdy_vvec;
00118     std::copy(bdy_vrange.begin(), bdy_vrange.end(), std::back_inserter(bdy_vvec));
00119     std::copy(new_ents.begin(), new_ents.end(), std::back_inserter(bdy_vvec));
00120     
00121       // now convert vertex indices into handles in-place, working from the back
00122     for (int i = 3*num_tris-1; i >= 0; i--) {
00123       assert(connecti[i] >= 0 && connecti[i] < (int)bdy_vvec.size());
00124       connect[i] = bdy_vvec[connecti[i]];
00125     }
00126     
00127       // put new tris into new entity range, then commit the mesh
00128     new_ents.insert(starth, starth+num_tris-1);
00129     me->commit_mesh(new_ents, COMPLETE_MESH);
00130   }
00131 }
00132 
00133 void CAMALTriAdvance::print_debug(ModelEnt *me, std::vector<double> &coords,
00134                                   moab::Range &bdy_vrange, std::vector<int> &group_sizes,
00135                                   std::vector<int> &bdy_ids)
00136 {
00137   std::cout << "Surface_bounadry_mesh: mesh_size = "
00138             << me->mesh_interval_size() << std::endl;
00139   
00140   for (int i = 0; i < (int) bdy_vrange.size(); i++) {
00141     std::cout << coords[3 * i] << "  " << coords[3 * i + 1]
00142               << "  " << coords[3 * i + 2] << std::endl;
00143   }
00144   
00145   std::cout << "bdy_vertex_size:" << bdy_vrange.size()
00146             << ", group_size:" << group_sizes.size() << std::endl;
00147   
00148   int index = 0;
00149   for (int i = 0; i < (int) group_sizes.size(); i++) {
00150     int g_size = group_sizes[i];
00151     std::cout << "boundary_order_group" << i + 1 << ", group_size="
00152               << g_size << std::endl;
00153     for (int j = 0; j < g_size; j++) {
00154           std::cout << bdy_ids[index + j] << " ";
00155     }
00156     std::cout << std::endl;
00157     index += g_size;
00158   }
00159   
00160   moab::ErrorCode rval;
00161   moab::EntityHandle outset;
00162   std::string outfile;
00163   std::stringstream ss;
00164   
00165   RefEntity* entity = reinterpret_cast<RefEntity*> (me->geom_handle());
00166   ss << "CAMALTri_boundary_surf";
00167   ss << entity->id();
00168   ss << "_";
00169 #ifdef PARALLEL
00170 #ifdef HAVE_PARALLEL_MOAB
00171   moab::ParallelComm* pcomm = moab::ParallelComm::get_pcomm(mk_core()->moab_instance(), 0);
00172   ss << "proc";
00173   ss << pcomm->proc_config().proc_rank();
00174 #endif
00175 #endif
00176   ss >> outfile;
00177   outfile += ".vtk";
00178   rval = mk_core()->moab_instance()->create_meshset(0, outset);
00179   MBERRCHK(rval, mk_core()->moab_instance());
00180   rval = mk_core()->moab_instance()->add_entities(outset, bdy_vrange);
00181   MBERRCHK(rval, mk_core()->moab_instance());
00182   rval = mk_core()->moab_instance()->write_mesh(outfile.c_str(), &outset, 1);
00183   MBERRCHK(rval, mk_core()->moab_instance());
00184   std::cout << outfile.c_str() << " is saved." << std::endl;
00185 }
00186 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines