MeshKit
1.0
|
00001 #include "meshkit/CAMALPaver.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 "CMLPaver.hpp" 00011 #include <vector> 00012 00013 namespace MeshKit 00014 { 00015 00016 moab::EntityType CAMALPaver::meshTps[] = {moab::MBVERTEX, moab::MBQUAD, moab::MBMAXTYPE}; 00017 00018 CAMALPaver::CAMALPaver(MKCore *mk_core, const MEntVector &me_vec) 00019 : MeshScheme(mk_core, me_vec) 00020 { 00021 } 00022 00023 CAMALPaver::~CAMALPaver() 00024 { 00025 } 00026 00027 void CAMALPaver::setup_this() 00028 { 00029 // need to constrain all edges to be even 00030 constrain_even(); 00031 00032 // then call setup_boundary, to set up edge meshers 00033 setup_boundary(); 00034 } 00035 00036 void CAMALPaver::execute_this() 00037 { 00038 00039 #ifdef HAVE_FBIGEOM 00040 if (mentSelection.empty()) 00041 { 00042 // create model ents from previous op 00043 create_model_ents_from_previous_ops(); 00044 // now look at the latest SizingFunction, and set it or each model ent 00045 int latestIndexSF = 0; // maybe we would need to set it right 00046 for (MEntSelection::iterator sit = mentSelection.begin(); 00047 sit != mentSelection.end(); sit++) { 00048 // make a me, for convenience 00049 ModelEnt *me = (*sit).first; 00050 me->sizing_function_index(latestIndexSF); // need to work on this one; how do we know? 00051 } 00052 // now, force setup of this node again, as we have added model entities to it 00053 setup_called(false); 00054 mk_core()->setup(false); 00055 // debug 00056 mk_core()->print_graph(); 00057 // it may not be enough, we may have to execute the previous ops, that were just 00058 // created during setup... not very clean code; 00059 mk_core()->execute_before((GraphNode *) this); 00060 } 00061 #endif 00062 for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) { 00063 // make a me, for convenience 00064 ModelEnt *me = (*sit).first; 00065 00066 if(me->get_meshed_state()>=COMPLETE_MESH) 00067 continue; 00068 // create a surface evaluator for this modelent, and a size evaluator 00069 CAMALSurfEval cse(me); 00070 00071 SizingFunction * sz = mk_core()->sizing_function(me->sizing_function_index()); 00072 if (!sz->variable()) 00073 sz = NULL; // no function variable 00074 00075 CAMALSizeEval mesize(me->mesh_interval_size(), sz); 00076 // make sure the size isn't negative 00077 // make sure the size isn't negative 00078 if (mesize.get_size() == -1) mesize.set_size(1.0); 00079 00080 // assemble bounding mesh 00081 std::vector<moab::EntityHandle> bdy; 00082 std::vector<int> group_sizes, bdy_ids; 00083 me->boundary(0, bdy, NULL, &group_sizes); 00084 00085 // convert the handles to integers for input to Paver 00086 moab::Range bdy_vrange; 00087 std::vector<double> coords; 00088 me->get_indexed_connect_coords(bdy, NULL, NULL, bdy_ids, coords, &bdy_vrange); 00089 00090 // now construct the CAMAL mesher, and pass it initial conditions 00091 CMLPaver paver(&cse, &mesize); 00092 bool success = paver.set_boundary_mesh(bdy_vrange.size(), &coords[0], group_sizes.size(), &group_sizes[0], &bdy_ids[0]); 00093 if (!success) 00094 ECERRCHK(MK_FAILURE, "Trouble setting boundary mesh."); 00095 00096 // ok, now generate the mesh 00097 int num_pts, num_quads; 00098 success = paver.generate_mesh(num_pts, num_quads); 00099 if (!success) 00100 ECERRCHK(MK_FAILURE, "Trouble generating quad mesh."); 00101 00102 moab::Range &new_ents = (*sit).second; 00103 moab::ErrorCode rval; 00104 00105 // resize the coords array, then get the coords of the new points 00106 assert(num_pts >= (int)bdy_vrange.size()); 00107 if (num_pts > (int)bdy_vrange.size()) { 00108 coords.resize(3*(num_pts-bdy_vrange.size())); 00109 int pts_returned = paver.get_points_buf(coords.size(), &coords[0], bdy_vrange.size()); 00110 if (pts_returned != num_pts-(int)bdy_vrange.size()) 00111 ECERRCHK(MK_FAILURE, "Number of new points returned from Paver doesn't agree with previous value output."); 00112 00113 // create the new vertices' entities on the face 00114 rval = mk_core()->moab_instance()->create_vertices(&coords[0], pts_returned, new_ents); 00115 MBERRCHK(rval, mk_core()->moab_instance()); 00116 } 00117 00118 // for quads, pre-allocate connectivity 00119 moab::ReadUtilIface *iface; 00120 rval = mk_core()-> moab_instance() -> query_interface(iface); 00121 MBERRCHK(rval, mk_core()->moab_instance()); 00122 00123 //create the quads, get a direct ptr to connectivity 00124 moab::EntityHandle starth, *connect; 00125 rval = iface->get_element_connect(num_quads, 4, moab::MBQUAD, 1, starth, connect); 00126 MBERRCHK(rval, mk_core()->moab_instance()); 00127 00128 // read connectivity directly into that array, as int's 00129 int *connecti = (int*) connect; 00130 int quads_returned = paver.get_quads_buf(4*num_quads, connecti); 00131 if (quads_returned != num_quads) 00132 ECERRCHK(MK_FAILURE, "Number of new quads returned from Paver doesn't agree with previous value output."); 00133 00134 // put vertex handles into an indexible array 00135 std::vector<moab::EntityHandle> bdy_vvec; 00136 std::copy(bdy_vrange.begin(), bdy_vrange.end(), std::back_inserter(bdy_vvec)); 00137 std::copy(new_ents.begin(), new_ents.end(), std::back_inserter(bdy_vvec)); 00138 00139 // now convert vertex indices into handles in-place, working from the back 00140 for (int i = 4*num_quads-1; i >= 0; i--) { 00141 assert(connecti[i] >= 0 && connecti[i] < (int)bdy_vvec.size()); 00142 connect[i] = bdy_vvec[connecti[i]]; 00143 } 00144 00145 // put new quads into new entity range, then commit the mesh 00146 new_ents.insert(starth, starth+num_quads-1); 00147 me->commit_mesh(new_ents, COMPLETE_MESH); 00148 } 00149 } 00150 00151 } // namespace MeshKit