MeshKit  1.0
CurveFacetMeshReader.cpp
Go to the documentation of this file.
00001 
00002 #include "meshkit/MKCore.hpp"
00003 #include "meshkit/CurveFacetMeshReader.hpp"
00004 #include "meshkit/ModelEnt.hpp"
00005 #include <iostream>
00006 #include <iGeom.h>
00007 #include <cmath>
00008 #include "MBTagConventions.hpp"
00009 #include "moab/Core.hpp"
00010 #include "moab/GeomTopoTool.hpp"
00011 
00012 namespace MeshKit
00013 {
00014 // Construction Function for CurveFacetMeshReader
00015 
00016 moab::EntityType CurveFacetMeshReader_types[] = {moab::MBVERTEX, moab::MBEDGE, moab::MBMAXTYPE};
00017 const moab::EntityType* CurveFacetMeshReader::output_types()
00018   { return CurveFacetMeshReader_types; }
00019 
00020 
00021 CurveFacetMeshReader::CurveFacetMeshReader(MKCore *mk_core, const MEntVector &ments)
00022   : MeshScheme(mk_core, ments)
00023 {
00024   mk = mk_core;
00025 }
00026 
00027 // Destructor Function for CurveFacetMeshReader
00028 CurveFacetMeshReader::~CurveFacetMeshReader()
00029 {
00030 }
00031 
00032 void CurveFacetMeshReader::setup_this()
00033 {
00034 
00035   //create a vertex mesher
00036   int meshop_dim = 0;
00037   MeshOp *vm = (MeshOp*) mk_core()->construct_meshop(meshop_dim);
00038 
00039   for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
00040     { 
00041       ModelEnt *me = mit->first;
00042       //do an initial check that the model entity is of the correct dimension
00043       if(me->dimension() != 1)
00044         {
00045           std::cout << "Found an entity that is of an incorrect dimension" << std::endl;
00046           mentSelection.erase(mit);
00047           continue;
00048         } 
00049 
00050      
00051       //make sure that its children (vertices) have been meshed
00052       MEntVector children; 
00053       me->get_adjacencies(0, children);
00054 
00055       for(MEntVector::iterator j = children.begin(); j != children.end(); j++)
00056         {
00057           ModelEnt *child = *j;
00058           if(child->is_meshops_list_empty()) vm->add_modelent(child);
00059         }
00060 
00061     }
00062 
00063   //make sure we're meshing the veritces first
00064   mk_core()->insert_node(vm, this, mk_core()->root_node());
00065   
00066   //set faceting parameters to default values if the are not already set. 
00067 
00068   //set the faceting_tolerance
00069   if(!facet_tol) facet_tol = 1e-4;
00070 
00071   //set the geom_reabs value
00072   if(!geom_res) geom_res = 1e-6;
00073 
00074 }
00075 
00076 void CurveFacetMeshReader::execute_this()
00077 {
00078  
00079   for(MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
00080     {
00081 
00082       ModelEnt *me = mit->first;
00083       
00084     
00085       //---------------------------------//
00086       //            FACET                //
00087       //---------------------------------//
00088       facet(me);
00089 
00090       // Not sure if this needs to be done yet
00091       //set_senses(me);
00092 
00093       me->set_meshed_state(COMPLETE_MESH);
00094     }
00095 }
00096 
00097 void CurveFacetMeshReader::facet(ModelEnt *curve)
00098 {
00099   iGeom::EntityHandle h = curve->geom_handle();
00100   iMesh::EntitySetHandle sh = IBSH(curve->mesh_handle());
00101 
00102   //get the facets for this curve/ref_edge
00103   std::vector<double> pnts;
00104   std::vector<int> conn;
00105   iGeom::Error errval = mk_core()->igeom_instance()->getFacets(h,facet_tol,pnts,conn);
00106   if (errval != iBase_SUCCESS) throw Error(MK_FAILURE, "Unable get the facets of the curve.");
00107   //create vector for keeping track of the vertices
00108   std::vector<iBase_EntityHandle> verts;
00109           
00110 
00111   // loop over the facet points
00112   for(unsigned int j = 0; j < pnts.size(); j+=3)
00113     {
00114       //create a new vertex handle 
00115       iBase_EntityHandle v;
00116       mk_core()->imesh_instance()->createVtx(pnts[j],pnts[j+1],pnts[j+2],v);
00117       verts.push_back(v);
00118     }
00119 
00120   //--------------------------------------------------//
00121   //check that the start and end vertices are within  //
00122   //the geom resabs of the facet verts                //
00123   //--------------------------------------------------//
00124      
00125   //get the curve vertices
00126   MEntVector end_verts;
00127   curve->get_adjacencies(0, end_verts);
00128   
00129   //get the entity handle of the start vertex
00130   std::vector<moab::EntityHandle> dum_handles;
00131   end_verts.front()->get_mesh(0,dum_handles);
00132   assert( 1 == dum_handles.size() );   // should only be one vertex in a vertex entity set
00133   iMesh::EntityHandle start_vtx = IBEH(dum_handles.front());
00134 
00135   dum_handles.clear();
00136 
00137   // do the same for the end vertex
00138   end_verts.back()->get_mesh(0, dum_handles);
00139   assert( 1 == dum_handles.size() );   // should only be one vertex in a vertex entity set
00140   iMesh::EntityHandle end_vtx = IBEH(dum_handles.front());
00141 
00142   dum_handles.clear();
00143 
00144   //std::cout << "Number of end_verts: " << end_verts.size() << std::endl;
00145 
00146   //special case for a point curve
00147   if(verts.size() < 2)
00148     {
00149       if( start_vtx != end_vtx || curve->measure() > geom_res)
00150         {
00151           std::cout << "Warning: no facetting for curve " << curve->id() << std::endl;
00152           return;
00153         }
00154       else
00155         {
00156           mk_core()->imesh_instance()->addEntToSet(start_vtx, sh);
00157           return;
00158         }
00159     }
00160 
00161   //check for a closed curve
00162   bool closed = (mvtx2mvtx_dist(verts.front(),verts.back()) < geom_res);
00163 
00164   // check that geometry and facetting agree that the curve is closed
00165   if ( closed != (start_vtx == end_vtx))
00166     std::cout << "Warning: topology and geometry inconsistent for possibly closed curve "
00167               << curve->id() << std::endl;
00168 
00169   //check the proximity of the front and end vertices to the start and end points, respectively
00170   if(vtx2vtx_dist(end_verts.front()->geom_handle(), verts.front()) > geom_res ||
00171      vtx2vtx_dist(end_verts.back()->geom_handle(), verts.back()) > geom_res)
00172     {
00173       //try reversing the points
00174       //std::reverse(verts.begin(), verts.end());
00175       
00176       //check again, if this time it fails, give a warning
00177       if(vtx2vtx_dist(end_verts.front()->geom_handle(), verts.front()) > geom_res ||
00178          vtx2vtx_dist(end_verts.back()->geom_handle(), verts.back()) > geom_res)
00179         {
00180           std::cerr << "Warning: vertices not at the ends of the curve " << curve->id() << std::endl;
00181         }
00182     }
00183 
00184 
00185   //now replace start and end facet points with the curve's child vert(s)
00186   
00187   //capture the beginning and end vertex handles for deletion 
00188   iMesh::EntityHandle front_vert_to_del = verts.front();
00189   iMesh::EntityHandle back_vert_to_del = verts.back();
00190   
00191   verts.front() = start_vtx;
00192   verts.back() = end_vtx;
00193   
00194   //delete the old vertices
00195   mk_core()->imesh_instance()->deleteEnt(front_vert_to_del);
00196   mk_core()->imesh_instance()->deleteEnt(back_vert_to_del);
00197     
00198   //now create the edges with the vertices we've chosen
00199 
00200   //vector to contain the edges
00201   std::vector<iBase_EntityHandle> edges;
00202   //loop over vertices and create edges
00203   for (unsigned int j = 0; j < verts.size()-1; j++)
00204     {
00205       //create edges
00206       iBase_EntityHandle e; 
00207       mk_core()->imesh_instance()->createEnt(iMesh_LINE_SEGMENT, &verts[j], 2, e);
00208       edges.push_back(e);
00209     }
00210 
00211   //remove duplicate vertex if the curve is closed
00212   if(closed && start_vtx == end_vtx) verts.pop_back(); 
00213 
00214   //add vertices and edges to the entity set
00215   mk_core()->imesh_instance()->addEntArrToSet(&verts[0], verts.size(), sh);
00216   mk_core()->imesh_instance()->addEntArrToSet(&edges[0], edges.size(), sh);
00217 }
00218 
00219 void CurveFacetMeshReader::set_senses(ModelEnt *curve)
00220 {
00221  
00222   //get the geom_handle for this curve
00223   iGeom::EntityHandle gh = curve->geom_handle();
00224 
00225   //get all surfaces adjacent to this curve
00226   MEntVector adj_surfs;
00227   curve->get_adjacencies(2, adj_surfs);
00228 
00229   MEntVector::iterator i; 
00230   std::vector<int> senses;
00231   std::vector<moab::EntityHandle> meshsets;
00232   for(i = adj_surfs.begin(); i !=adj_surfs.end(); i++)
00233     {
00234       int sense;
00235       //create a model ent for the current adjacent surface
00236       ModelEnt *adj_ent = *i;
00237       //get the sense of the curve wrt the surface
00238       mk_core()->igeom_instance()->getEgFcSense(gh, adj_ent->geom_handle(),sense);
00239       //std::cout << "Sense: " << sense << std::endl;
00240       // add the sense and the entityhandle to the appropriate vectors
00241       senses.push_back(sense);
00242       meshsets.push_back(adj_ent->mesh_handle());
00243     }
00244 
00245   // use moab to set the senses of the entities
00246   // using moab to do this allows the use of variable length tags
00247   // which are not currently available in iGeom 
00248   moab::GeomTopoTool gt(mk_core()->moab_instance());
00249   moab::ErrorCode rval = gt.set_senses(curve->mesh_handle(),meshsets,senses);
00250   if(rval != moab::MB_SUCCESS) std::cout << "Error setting curve senses!" << std::endl;
00251 
00252 }
00253 
00254 // should I be using an iMesh entity to do this comparison??
00255 double CurveFacetMeshReader::vtx2vtx_dist(iGeom::EntityHandle vtx1, iMesh::EntityHandle vtx2)
00256 {
00257 
00258   double x1,y1,z1;
00259   double x2,y2,z2;
00260 
00261   mk_core()->igeom_instance()->getVtxCoord(vtx1, x1, y1, z1);
00262   mk_core()->imesh_instance()->getVtxCoord(vtx2, x2, y2, z2);
00263 
00264   double dist = pow((x1-x2),2) + pow((y1-y2),2) + pow((z1-z2),2);
00265 
00266   return sqrt(dist);
00267 }
00268 
00269 double CurveFacetMeshReader::mvtx2mvtx_dist(iMesh::EntityHandle vtx1, iMesh::EntityHandle vtx2)
00270 {
00271 
00272   double x1,y1,z1;
00273   double x2,y2,z2;
00274 
00275   mk_core()->imesh_instance()->getVtxCoord(vtx1, x1, y1, z1);
00276   mk_core()->imesh_instance()->getVtxCoord(vtx2, x2, y2, z2);
00277 
00278   double dist = pow((x1-x2),2) + pow((y1-y2),2) + pow((z1-z2),2);
00279 
00280   return sqrt(dist);
00281 }
00282 
00283 void CurveFacetMeshReader::set_mesh_params(double faceting_tolerance, double geom_resabs)
00284 {
00285   //Assign the faceting values if they are passed into the function
00286   if(faceting_tolerance) facet_tol = faceting_tolerance; 
00287   if(geom_resabs) geom_res = geom_resabs; 
00288 }
00289 
00290 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines