MeshKit
1.0
|
00001 #include "meshkit/OneToOneSwept.hpp" 00002 #include "meshkit/iMesh.hpp" 00003 #include "meshkit/EdgeMesher.hpp" 00004 #include "meshkit/TFIMapping.hpp" 00005 #ifdef HAVE_MESQUITE 00006 #include "meshkit/MeshImprove.hpp" 00007 #endif 00008 #include <iostream> 00009 #include <algorithm> 00010 #include <math.h> 00011 #include <map> 00012 #ifdef HAVE_ARMADILLO 00013 #include "SurfProHarmonicMap.hpp" 00014 #endif 00015 00016 namespace MeshKit { 00017 00018 //---------------------------------------------------------------------------// 00019 //Entity Type initialization for OneToOneSwept meshing 00020 moab::EntityType OneToOneSwept_tps[] = { moab::MBVERTEX, moab::MBQUAD, moab::MBHEX, moab::MBMAXTYPE }; 00021 const moab::EntityType* OneToOneSwept::output_types() 00022 { 00023 return OneToOneSwept_tps; 00024 } 00025 00026 //---------------------------------------------------------------------------// 00027 // construction function for OneToOneSwept class 00028 OneToOneSwept::OneToOneSwept(MKCore *mk_core, const MEntVector &me_vec) : 00029 MeshScheme(mk_core, me_vec) 00030 { 00031 } 00032 00033 //---------------------------------------------------------------------------// 00034 // deconstruction function for OneToOneSwept class 00035 OneToOneSwept::~OneToOneSwept() 00036 { 00037 } 00038 00039 // these have to be called before setup, otherwise we don't know source and target 00040 //---------------------------------------------------------------------------// 00041 // set the source surface function 00042 void OneToOneSwept::SetSourceSurface(int index) 00043 { 00044 index_src = index; 00045 } 00046 00047 //---------------------------------------------------------------------------// 00048 // set the target surface function 00049 void OneToOneSwept::SetTargetSurface(int index) 00050 { 00051 index_tar = index; 00052 } 00053 00054 //---------------------------------------------------------------------------// 00055 // setup function: define the size between the different layers 00056 void OneToOneSwept::setup_this() 00057 { 00058 //compute the number of intervals for the associated ModelEnts, from the size set on them 00059 //the sizing function they point to, or a default sizing function 00060 00061 if (mentSelection.empty()) 00062 return; 00063 ModelEnt * firstME = (mentSelection.begin())->first; 00064 igeom_inst = firstME->igeom_instance(); 00065 /*moab::Interface **/ 00066 mb = mk_core()->moab_instance();// we should get it also from 00067 // model entity, if we have multiple moab instances running around 00068 00069 const char *tag = "GLOBAL_ID"; 00070 iGeom::Error g_err = igeom_inst->getTagHandle(tag, geom_id_tag); 00071 IBERRCHK(g_err, "Trouble get the geom_id_tag for 'GLOBAL_ID'."); 00072 00073 for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++) 00074 { 00075 ModelEnt *me = mit->first; 00076 //first check to see whether entity is meshed 00077 if (me->get_meshed_state() >= COMPLETE_MESH || me->mesh_intervals() > 0) 00078 continue; 00079 00080 SizingFunction *sf = mk_core()->sizing_function(me->sizing_function_index()); 00081 if (!sf && me -> mesh_intervals() < 0 && me -> interval_firmness() == DEFAULT && mk_core()->sizing_function(0)) 00082 sf = mk_core()->sizing_function(0); 00083 00084 if (!sf && me -> mesh_intervals() < 0 && me -> interval_firmness() == DEFAULT) 00085 { 00086 //no sizing set, just assume default #intervals as 4 00087 me->mesh_intervals(4); 00088 me->interval_firmness(DEFAULT); 00089 } 00090 else 00091 { 00092 //check # intervals first, then size, and just choose for now 00093 if (sf->intervals() > 0) 00094 { 00095 if (me->constrain_even() && sf->intervals() % 2) 00096 me -> mesh_intervals(sf->intervals() + 1); 00097 else 00098 me -> mesh_intervals(sf->intervals()); 00099 me -> interval_firmness(HARD); 00100 } 00101 else if (sf->size() > 0) 00102 { 00103 int intervals = me->measure() / sf->size(); 00104 if (!intervals) 00105 intervals++; 00106 if (me->constrain_even() && intervals % 2) 00107 intervals++; 00108 me->mesh_intervals(intervals); 00109 me->interval_firmness(SOFT); 00110 } 00111 else 00112 throw Error(MK_INCOMPLETE_MESH_SPECIFICATION, "Sizing function for edge had neither positive size nor positive intervals."); 00113 } 00114 00115 int numIntervals = me->mesh_intervals(); 00116 if (!sf || sf->intervals() != numIntervals) 00117 { 00118 sf = new SizingFunction(mk_core(), numIntervals, -1.); 00119 } 00120 00121 std::vector<iBase_EntityHandle> gFaces; 00122 g_err = igeom_inst->getEntAdj(me->geom_handle(), iBase_FACE, gFaces); 00123 IBERRCHK(g_err, "Trouble get the geometrical adjacent to the solid"); 00124 //select the source surface and target surface 00125 // these were setup from the beginning, it could be different for each volume!!!! 00126 // we need a better way to decide/select target and source surfaces 00127 // this also assumes that the sourceSurface is fixed between setup and execute 00128 sourceSurface = gFaces[index_src]; 00129 targetSurface = gFaces[index_tar]; 00130 MEntVector linkingSurfaces; 00131 int index_id_src, index_id_tar; 00132 iGeom::Error g_err = igeom_inst->getIntData(sourceSurface, geom_id_tag, index_id_src); 00133 IBERRCHK(g_err, "Trouble get the int data for source surface."); 00134 g_err = igeom_inst->getIntData(targetSurface, geom_id_tag, index_id_tar); 00135 IBERRCHK(g_err, "Trouble get the int data for target surface."); 00136 std::cout << " source Global ID: " << index_id_src << " target Global ID: " << index_id_tar << "\n"; 00137 00138 for (unsigned int i = 0; i < gFaces.size(); i++) 00139 { 00140 int gid; 00141 g_err = igeom_inst->getIntData(gFaces[i], geom_id_tag, gid); 00142 IBERRCHK(g_err, "Trouble get the int data for source surface."); 00143 std::vector<iBase_EntityHandle> gEdges; 00144 g_err = igeom_inst->getEntAdj(gFaces[i], iBase_EDGE, gEdges); 00145 IBERRCHK(g_err, "Trouble get the geometrical edges adjacent to the surface"); 00146 std::cout << " face index " << i << " with ID: " << gid << " has " << gEdges.size() << " edges\n"; 00147 } 00148 MEntVector surfs; 00149 me->get_adjacencies(2, surfs);// all surfaces adjacent to the volume 00150 for (unsigned int i = 0; i < surfs.size(); i++) 00151 { 00152 int index_id_link; 00153 g_err = igeom_inst->getIntData(surfs[i]->geom_handle(), geom_id_tag, index_id_link); 00154 IBERRCHK(g_err, "Trouble get the int data for linking surface."); 00155 if ((index_id_link != index_id_src) && (index_id_link != index_id_tar)) 00156 { 00157 MEntVector edges; 00158 surfs[i]->get_adjacencies(1, edges);// all surfaces adjacent to the volume 00159 std::cout << " linking surface with global ID: " << index_id_link << " has " << edges.size() << " edges\n"; 00160 assert((int)edges.size()==4); 00161 00162 linkingSurfaces.push_back(surfs[i]); 00163 } 00164 } 00165 // now create a TFI mapping 00166 // this will trigger edge meshers for linking edges, and for target surface edges 00167 TFIMapping *tm = (TFIMapping*) mk_core()->construct_meshop("TFIMapping", linkingSurfaces); 00168 for (unsigned int i = 0; i < linkingSurfaces.size(); i++) 00169 { 00170 linkingSurfaces[i]->sizing_function_index(sf->core_index()); 00171 } 00172 mk_core()->insert_node(tm, (MeshOp*) this); 00173 } 00174 00175 mk_core()->print_graph("AfterOneSetup.eps"); 00176 00177 } 00178 00179 //---------------------------------------------------------------------------// 00180 // execute function: generate the all-hex mesh through sweeping from source 00181 // surface to target surface 00182 void OneToOneSwept::execute_this() 00183 { 00184 00185 for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++) 00186 { 00187 ModelEnt *me = mit -> first; 00188 if (me->get_meshed_state() >= COMPLETE_MESH) 00189 continue; 00190 numLayers = me->mesh_intervals();// maybe it will be decided later? 00191 //case 1: if numLayers = 1, then it is not necessary to create any vertices for linking surface, All the vertices have been created by source surface and target surface 00192 // bLayers will contain nodes in order, layer by layer 00193 std::vector<moab::EntityHandle> bLayers; 00194 BuildLateralBoundaryLayers(me, bLayers); 00195 //int sizeBLayer = bLayers.size()/(numLayers+1); 00196 volume = me->geom_handle(); 00197 00198 TargetSurfProjection(bLayers);// the top layer is the last list of nodes 00199 00200 //get the volume mesh entity set 00201 iRel::Error r_err = mk_core()->irel_pair(me->iRelPairIndex())->getEntSetRelation(me->geom_handle(), 0, volumeSet); 00202 IBERRCHK(r_err, 00203 "Trouble get the volume mesh entity set from the geometrical volume."); 00204 00205 vector<vector<Vertex> > linkVertexList(numLayers - 1, vector<Vertex> (NodeList.size())); 00206 // this will compute the interior points distribution, in each layer 00207 ProjectInteriorLayers(bLayers, linkVertexList); 00208 00209 CreateElements(linkVertexList); 00210 00211 mk_core()->save_mesh("BeforeVolumeImprove.h5m"); 00212 #if HAVE_MESQUITE 00213 MeshImprove meshImpr(mk_core()); 00214 meshImpr.VolumeMeshImprove(volumeSet, iBase_REGION); 00215 #endif 00216 me->commit_mesh(mit->second, COMPLETE_MESH); 00217 } 00218 00219 } 00220 00221 // will also initialize NodeList, which comprises all <Vertex> data on the source sf 00222 void OneToOneSwept::BuildLateralBoundaryLayers(ModelEnt * me, std::vector<moab::EntityHandle> & layers) 00223 { 00224 // start with the source surface 00225 //ModelEnt * sourceME = ModelEnt() 00226 iGeom * igeom_inst = me->igeom_instance(); 00227 iGeom::Error g_err; 00228 moab::ErrorCode rval = moab::MB_SUCCESS; 00229 int index_id_src, index_id_tar; 00230 00231 g_err = igeom_inst->getIntData(sourceSurface, geom_id_tag, index_id_src); 00232 IBERRCHK(g_err, "Trouble get the int data for source surface."); 00233 g_err = igeom_inst->getIntData(targetSurface, geom_id_tag, index_id_tar); 00234 IBERRCHK(g_err, "Trouble get the int data for target surface."); 00235 00236 iBase_TagHandle taghandle = 0; 00237 iMesh::Error m_err = mk_core()->imesh_instance()->createTag("source", 1, iBase_INTEGER, taghandle); 00238 IBERRCHK(m_err, "Trouble create the tag handle in the mesh instance."); 00239 00240 MEntVector surfs; 00241 me->get_adjacencies(2, surfs);// all surfaces adjacent to the volume 00242 00243 moab::Range quadsOnLateralSurfaces; 00244 00245 ModelEnt * sME = NULL; 00246 me->get_adjacencies(2, surfs);// all surfaces adjacent to the volume 00247 for (unsigned int i = 0; i < surfs.size(); i++) 00248 { 00249 int index_id_link; 00250 g_err = igeom_inst->getIntData(surfs[i]->geom_handle(), geom_id_tag, index_id_link); 00251 IBERRCHK(g_err, "Trouble get the int data for linking surface."); 00252 if ((index_id_link != index_id_src) && (index_id_link != index_id_tar)) 00253 { 00254 moab::EntityHandle surfSet = surfs[i]->mesh_handle(); 00255 // get all nodes from the surface, including the boundary nodes! 00256 rval = mb->get_entities_by_type(surfSet, moab::MBQUAD, quadsOnLateralSurfaces); 00257 MBERRCHK(rval, mb); 00258 } 00259 if (index_id_link == index_id_src) 00260 sME = surfs[i]; 00261 } 00262 // the lateral layers will be build from these nodes 00263 moab::Range nodesOnLateralSurfaces; 00264 rval = mb->get_connectivity(quadsOnLateralSurfaces, nodesOnLateralSurfaces); 00265 MBERRCHK(rval, mb); 00266 00267 // construct NodeList, nodes on source surface 00268 // list of node on boundary of source 00269 std::vector<moab::EntityHandle> bdyNodes; 00270 std::vector<int> group_sizes; 00271 sME-> boundary(0, bdyNodes, NULL, &group_sizes);// we do not need the senses, yet 00272 00273 numLoops = (int) group_sizes.size(); 00274 sizeBLayer = (int) bdyNodes.size(); 00275 00276 std::cout << "Execute one to one swept ....\n"; 00277 std::cout << "Nodes on boundary: " << bdyNodes.size() << "\n"; 00278 std::cout << "Number of loops: " << group_sizes.size() << "\n"; 00279 std::cout << "Nodes on lateral surfaces: " << nodesOnLateralSurfaces.size() << "\n"; 00280 std::cout << "Quads on lateral surfaces: " << quadsOnLateralSurfaces.size() << "\n"; 00281 00282 if (bdyNodes.size() * (numLayers + 1) != nodesOnLateralSurfaces.size()) 00283 { 00284 std::cout << "Major problem: number of total nodes on boundary: " << nodesOnLateralSurfaces.size() << "\n"; 00285 std::cout << " number of nodes in one layer on the boundary:" << bdyNodes.size() << "\n"; 00286 std::cout << " we expect bdyNodes.size()*(numLayers+1) == nodesOnLateralSurfaces.size()" << bdyNodes.size() * (numLayers + 1) 00287 << " != " << nodesOnLateralSurfaces.size() << "\n"; 00288 } 00289 // start from this list of boundary nodes (on the source) 00290 00291 // get all nodes from the source surface 00292 moab::Range sourceNodes; 00293 moab::EntityHandle surfSet = sME->mesh_handle(); 00294 // get all nodes from the surface, including the boundary nodes! 00295 moab::Range quads; 00296 rval = mb->get_entities_by_type(surfSet, moab::MBQUAD, quads); 00297 MBERRCHK(rval, mb); 00298 // we do not really care about corners, only about boundary 00299 rval = mb->get_connectivity(quads, sourceNodes); 00300 MBERRCHK(rval, mb); 00301 00302 NodeList.resize(sourceNodes.size()); 00303 int i = 0; 00304 for (moab::Range::iterator it = sourceNodes.begin(); it != sourceNodes.end(); it++, i++) 00305 { 00306 moab::EntityHandle node = *it; 00307 NodeList[i].gVertexHandle = (iBase_EntityHandle) sourceNodes[i]; 00308 NodeList[i].index = i; 00309 00310 rval = mb->get_coords(&node, 1, &(NodeList[i].xyz[0])); 00311 MBERRCHK(rval, mb); 00312 NodeList[i].onBoundary = false; 00313 00314 //set the int data to the entity 00315 m_err = mk_core()->imesh_instance()->setIntData(NodeList[i].gVertexHandle, taghandle, NodeList[i].index); 00316 IBERRCHK(m_err, "Trouble set the int value for nodes in the mesh instance."); 00317 } 00318 // now loop over the boundary and mark the boundary nodes as such 00319 for (unsigned int j = 0; j < bdyNodes.size(); j++) 00320 { 00321 iBase_EntityHandle node = (iBase_EntityHandle) bdyNodes[j]; 00322 int index; 00323 m_err = mk_core()->imesh_instance()->getIntData(node, taghandle, index); 00324 IBERRCHK(m_err, "Trouble set the int value for nodes in the mesh instance."); 00325 NodeList[index].onBoundary = true;// some could be corners, but it is not that important 00326 } 00327 00328 // process faces on the source 00329 FaceList.resize(quads.size());// quads is actually a range.... 00330 00331 for (unsigned int i = 0; i < quads.size(); i++) 00332 { 00333 iBase_EntityHandle currentFace = (iBase_EntityHandle) quads[i];// 00334 FaceList[i].gFaceHandle = currentFace; 00335 //FaceList[i].index = i; 00336 00337 //get the nodes on the face elements 00338 std::vector<iBase_EntityHandle> Nodes; 00339 m_err = mk_core()->imesh_instance()->getEntAdj(currentFace, iBase_VERTEX, Nodes); 00340 IBERRCHK(m_err, "Trouble get the adjacent nodes from mesh face entity."); 00341 00342 FaceList[i].connect.resize(Nodes.size()); 00343 for (unsigned int j = 0; j < Nodes.size(); j++) 00344 { 00345 int tmpIndex = -1; 00346 m_err = mk_core()->imesh_instance()->getIntData(Nodes[j], taghandle, tmpIndex); 00347 IBERRCHK(m_err, "Trouble get the int data from node handle."); 00348 00349 FaceList[i].connect[j] = &NodeList[tmpIndex];// why not store tmp Index directly? 00350 } 00351 m_err = mk_core()->imesh_instance()->setIntData(currentFace, taghandle, i); 00352 IBERRCHK(m_err, "Trouble set the int data for quads on the source surface."); 00353 00354 /* m_err = mk_core()->imesh_instance()->setIntData(currentFace, taghandle, i); 00355 IBERRCHK(m_err, "Trouble set the int data for quad mesh on the source surface.");*/ 00356 } 00357 00358 std::copy(bdyNodes.begin(), bdyNodes.end(), std::back_inserter(layers)); 00359 00360 int deft = -1; // not set; 0 mean layer 1; maybe; we should maybe not use it? 00361 rval = mb->tag_get_handle("mark", 1, moab::MB_TYPE_INTEGER, markTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &deft); 00362 MBERRCHK(rval, mb); 00363 00364 int val = 0;// layer level ; we will use this to verify the layers of boundary nodes 00365 for (i = 0; i < sizeBLayer; i++) 00366 { 00367 rval = mb->tag_set_data(markTag, &(bdyNodes[i]), 1, &val); 00368 MBERRCHK(rval, mb); 00369 } 00370 00371 for (int layer = 1; layer <= numLayers; layer++) // layer with index 0 is at the bottom/source 00372 { 00373 int start_index = 0;//this is the start index in group 00374 for (unsigned int k = 0; k < group_sizes.size(); k++) 00375 { 00376 // first group starts at index 0, the rest are accumulated 00377 if (k > 0) 00378 start_index += group_sizes[k - 1]; 00379 moab::EntityHandle node1 = layers[(layer - 1) * sizeBLayer + start_index]; 00380 moab::EntityHandle node2 = layers[(layer - 1) * sizeBLayer + start_index + 1]; 00381 moab::EntityHandle node3, node4; 00382 /* 00383 * node3 -- node4 --> (next layer) 00384 * | | 00385 * node1 -- node2 ---> (guide layer) 00386 */ 00387 rval = NodeAbove(node1, node2, quadsOnLateralSurfaces, node3, node4); 00388 MBERRCHK(rval, mb); 00389 // check the mark tag! 00390 rval = mb->tag_get_data(markTag, &node3, 1, &val); 00391 MBERRCHK(rval, mb); 00392 if (val != -1) 00393 MBERRCHK(moab::MB_FAILURE, mb); 00394 rval = mb->tag_get_data(markTag, &node4, 1, &val); 00395 MBERRCHK(rval, mb); 00396 if (val != -1) 00397 MBERRCHK(moab::MB_FAILURE, mb); 00398 layers.push_back(node3);//start of a new layer 00399 layers.push_back(node4);//node above node 2 00400 rval = mb->tag_set_data(markTag, &node3, 1, &layer);// layer level>0 00401 MBERRCHK(rval, mb); 00402 rval = mb->tag_set_data(markTag, &node4, 1, &layer);// layer level>0 00403 MBERRCHK(rval, mb); 00404 // we have at least 2 nodes in every loop, otherwise we are in deep trouble 00405 //start 00406 int currentIndex = start_index + 2; 00407 while (currentIndex < start_index + group_sizes[k]) 00408 { 00409 node1 = node2; 00410 node2 = layers[(layer - 1) * sizeBLayer + currentIndex];// actually, the previous layer 00411 node3 = node4; 00412 rval = FourthNodeInQuad(node1, node2, node3, quadsOnLateralSurfaces, node4); 00413 MBERRCHK(rval, mb); 00414 rval = mb->tag_get_data(markTag, &node4, 1, &val); 00415 MBERRCHK(rval, mb); 00416 if (val != -1) 00417 MBERRCHK(moab::MB_FAILURE, mb); 00418 layers.push_back(node4); 00419 rval = mb->tag_set_data(markTag, &node4, 1, &layer);// layer level>0 00420 MBERRCHK(rval, mb); 00421 currentIndex++; 00422 } 00423 }// end group 00424 }// end layer 00425 for (int lay = 0; lay <= numLayers; lay++) 00426 { 00427 std::cout << "layer : " << lay << "\n"; 00428 for (int i = 0; i < sizeBLayer; i++) 00429 { 00430 std::cout << " " << layers[lay * sizeBLayer + i]; 00431 if (i % 10 == 9) 00432 std::cout << "\n"; 00433 } 00434 std::cout << "\n"; 00435 } 00436 } 00437 00438 moab::ErrorCode OneToOneSwept::NodeAbove(moab::EntityHandle node1, moab::EntityHandle node2, moab::Range & latQuads, 00439 moab::EntityHandle & node3, moab::EntityHandle & node4) 00440 { 00441 // look for node above node 1 in an unused quad 00442 moab::EntityHandle nds[2] = { node1, node2 }; 00443 // find all quads connected to these 2 nodes; find one in the range that is not used 00444 moab::Range adj_entities; 00445 moab::ErrorCode rval = mb->get_adjacencies(nds, 2, 2, false, adj_entities); 00446 MBERRCHK(rval, mb); 00447 std::vector<int> tagVals(adj_entities.size()); 00448 rval = mb->tag_get_data(markTag, adj_entities, &(tagVals[0])); 00449 MBERRCHK(rval, mb); 00450 // default is -1 00451 moab::EntityHandle nextQuad = 0;// null 00452 for (unsigned int i = 0; i < tagVals.size(); i++) 00453 { 00454 if ((tagVals[i] == -1) && (latQuads.find(adj_entities[i]) != latQuads.end())) 00455 { 00456 nextQuad = adj_entities[i]; 00457 break; 00458 } 00459 } 00460 if (0 == nextQuad) 00461 MBERRCHK(moab::MB_FAILURE, mb); 00462 00463 // decide on which side are nodes 1 and 2, then get the opposite side 00464 const moab::EntityHandle * conn4 = NULL; 00465 int nnodes; 00466 rval = mb->get_connectivity(nextQuad, conn4, nnodes); 00467 MBERRCHK(rval, mb); 00468 int index1 = -1; 00469 for (index1 = 0; index1 < 4; index1++) 00470 if (node1 == conn4[index1]) 00471 break; 00472 if (4 == index1) 00473 MBERRCHK(moab::MB_FAILURE, mb); 00474 if (node2 == conn4[(index1 + 1) % 4]) // quad is oriented node1, node2, node4, node3 00475 { 00476 node3 = conn4[(index1 + 3) % 4]; 00477 } 00478 else if (node2 == conn4[(index1 +3) % 4]) // quad is oriented node2, node1, node3, node4 00479 { 00480 node3 = conn4[(index1 + 1) % 4]; 00481 } 00482 else 00483 MBERRCHK(moab::MB_FAILURE, mb);// something is really wrong 00484 node4 = conn4[(index1 + 2) % 4]; 00485 00486 // mark the quad to not use it again 00487 int val = 0; // maybe it should be the layer 00488 rval = mb->tag_set_data(markTag, &nextQuad, 1, &val); 00489 MBERRCHK(rval, mb); 00490 00491 return moab::MB_SUCCESS; 00492 } 00493 00494 moab::ErrorCode OneToOneSwept::FourthNodeInQuad(moab::EntityHandle node1, moab::EntityHandle node2, moab::EntityHandle node3, 00495 moab::Range & latQuads, moab::EntityHandle & node4) 00496 { 00497 // look for the fourth node 00498 moab::EntityHandle nds[3] = { node1, node2, node3 }; 00499 // find all quads connected to these 3 nodes; find one in the range that is not used 00500 moab::Range adj_entities; 00501 moab::ErrorCode rval = mb->get_adjacencies(nds, 3, 2, false, adj_entities); 00502 MBERRCHK(rval, mb); 00503 00504 std::vector<int> tagVals(adj_entities.size()); 00505 rval = mb->tag_get_data(markTag, adj_entities, &(tagVals[0])); 00506 00507 MBERRCHK(rval, mb); 00508 // default is -1 00509 moab::EntityHandle nextQuad = 0;// null 00510 for (unsigned int i = 0; i < tagVals.size(); i++) 00511 { 00512 if ((tagVals[i] == -1) && (latQuads.find(adj_entities[i]) != latQuads.end())) 00513 { 00514 nextQuad = adj_entities[i]; 00515 break; 00516 } 00517 } 00518 if (0 == nextQuad) 00519 MBERRCHK(moab::MB_FAILURE, mb); 00520 00521 // decide on which side are nodes 1 and 2, then get the opposite side 00522 const moab::EntityHandle * conn4 = NULL; 00523 int nnodes; 00524 rval = mb->get_connectivity(nextQuad, conn4, nnodes); 00525 MBERRCHK(rval, mb); 00526 00527 node4 = 0; 00528 for (int index = 0; index < 4; index++) 00529 { 00530 moab::EntityHandle c4 = conn4[index]; 00531 if (node1 != c4 && node2 != c4 && node3 != c4) 00532 { 00533 node4 = c4; 00534 break; 00535 } 00536 } 00537 if (0 == node4) 00538 MBERRCHK(moab::MB_FAILURE, mb); 00539 00540 // mark the quad to not use it again 00541 int val = 0; // maybe it should be the layer 00542 rval = mb->tag_set_data(markTag, &nextQuad, 1, &val); 00543 MBERRCHK(rval, mb); 00544 00545 return moab::MB_SUCCESS; 00546 00547 } 00548 #ifdef HAVE_ARMADILLO 00549 void OneToOneSwept::SurfMeshHarmonic(iBase_EntityHandle vol, vector<iBase_EntityHandle> &newNodehandle) 00550 { 00551 00552 SurfProHarmonicMap surf_pro(mk_core(), sourceSurface, targetSurface, vol); 00553 surf_pro.match(); 00554 surf_pro.setMeshData(NodeList, TVertexList, FaceList); 00555 surf_pro.projection(); 00556 surf_pro.getMeshData(TVertexList); 00557 00558 iMesh::Error m_err; 00559 for (unsigned int i = 0; i < TVertexList.size(); i++){ 00560 if (TVertexList[i].onBoundary) 00561 continue; 00562 m_err = mk_core()->imesh_instance()->createVtx(TVertexList[i].xyz[0], TVertexList[i].xyz[1], TVertexList[i].xyz[2], TVertexList[i].gVertexHandle); 00563 IBERRCHK(m_err, "Trouble create a mesh nodes on the target surface!"); 00564 newNodehandle.push_back(TVertexList[i].gVertexHandle); 00565 } 00566 } 00567 #endif 00568 00569 bool OneToOneSwept::isConcave() 00570 { 00571 //use the quad mesh on the source surface to determine whether it is concave or not 00572 iBase_TagHandle taghandle = 0; 00573 iMesh::Error m_err = mk_core()->imesh_instance()->getTagHandle("source", taghandle); 00574 IBERRCHK(m_err, "Trouble get the tag handle in the mesh instance."); 00575 for (std::vector<Vertex>::iterator it = NodeList.begin(); it != NodeList.end(); it++){ 00576 if (!it->onBoundary) 00577 continue; 00578 //get the adjacent quads around a vertex 00579 double angle = 0.0; 00580 std::vector<iBase_EntityHandle> adj_quads; 00581 m_err = mk_core()->imesh_instance()->getEntAdj(it->gVertexHandle, iBase_FACE, adj_quads); 00582 IBERRCHK(m_err, "Trouble get the adjacent quads around a vertex!"); 00583 for (unsigned int i = 0; i < adj_quads.size(); i++){ 00584 int face_index = -1; 00585 m_err = mk_core()->imesh_instance()->getIntData(adj_quads[i], taghandle, face_index); 00586 if (m_err) continue;//this is a quad entity from the linking surface 00587 int node_index, pre_index, next_index; 00588 for (int j = 0; j < 4; j++) 00589 if (FaceList[face_index].connect[j]->index == it->index){ 00590 node_index = j; 00591 break; 00592 } 00593 Vector3D v1(0.0), v2(0.0); 00594 pre_index = (FaceList[face_index].getVertex((node_index+4-1)%4))->index; 00595 next_index = (FaceList[face_index].getVertex((node_index+1)%4))->index; 00596 v1 = NodeList[pre_index].xyz - it->xyz; 00597 v2 = NodeList[next_index].xyz - it->xyz; 00598 double dotproduct = v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]; 00599 angle += acos(dotproduct/sqrt((v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2])*(v2[0]*v2[0]+v2[1]*v2[1]+v2[2]*v2[2]))); 00600 } 00601 if (angle > PI){ 00602 return true; 00603 } 00604 } 00605 00606 return false; 00607 } 00608 00609 00610 //****************************************************************************// 00611 // function : TargetSurfProjection 00612 // Author : Shengyong Cai 00613 // Date : Feb 15, 2011 00614 // Description: map the mesh on the source surface to the target surface 00615 //***************************************************************************// 00616 int OneToOneSwept::TargetSurfProjection(std::vector<moab::EntityHandle> & bLayers) 00617 { 00618 iMesh::Error m_err; 00619 //iGeom::Error g_err; 00620 iRel::Error r_err; 00621 moab::ErrorCode rval; 00622 00623 MEntVector surfs; 00624 mk_core()->get_entities_by_dimension(2, surfs); 00625 ModelEnt *target_surf = 0; 00626 for (unsigned int i = 0; i < surfs.size(); i++) 00627 { 00628 00629 if (surfs[i]->geom_handle() == targetSurface) 00630 { 00631 target_surf = surfs[i]; 00632 break; 00633 } 00634 } 00635 int irelIndx = target_surf->iRelPairIndex(); 00636 iGeom * igeom_inst = mk_core()->igeom_instance(target_surf->iGeomIndex()); 00637 // we could have got it from model tag, too 00638 // something along this: 00639 // err = igeom_instance(geom_index)->getData(*eit, iGeomModelTags[geom_index], &this_me); 00640 00641 TVertexList.resize(NodeList.size()); 00642 // make everything not on boundary, first; the true boundary flag will be set later 00643 for (unsigned int i = 0; i < TVertexList.size(); i++) 00644 { 00645 TVertexList[i].onBoundary = false; 00646 } 00647 00648 iBase_TagHandle taghandle_tar = 0; 00649 m_err = mk_core()->imesh_instance()->createTag("TargetMesh", 1, iBase_INTEGER, taghandle_tar); 00650 IBERRCHK(m_err, "Trouble create the taghandle for the target mesh."); 00651 00652 iBase_TagHandle taghandle = 0; 00653 m_err = mk_core()->imesh_instance()->getTagHandle("source", taghandle); 00654 IBERRCHK(m_err, "Trouble get tag handle of the source surface."); 00655 00656 // we know the first layer of nodes (on the source) 00657 // 00658 // the first 0 , 1, ..., sizeBlayer - 1 are on bottom, were bdyNodes before 00659 // 00660 for (int i = 0; i < sizeBLayer; i++) 00661 { 00662 iBase_EntityHandle sBNode = (iBase_EntityHandle) bLayers[i]; 00663 int index; 00664 m_err = mk_core()->imesh_instance()->getIntData(sBNode, taghandle, index); 00665 IBERRCHK(m_err, "Trouble get the int value for nodes in the mesh instance."); 00666 TVertexList[index].onBoundary = true;// some could be corners, but it is not that important 00667 moab::EntityHandle topNode = bLayers[numLayers * sizeBLayer + i];// node right above 00668 iBase_EntityHandle tBNode = (iBase_EntityHandle) topNode;// it is right above node i in loop 00669 TVertexList[index].gVertexHandle = tBNode; 00670 // now get the coordinates of the tBNode (node on boundary of target) 00671 rval = mb->get_coords(&topNode, 1, &(TVertexList[index].xyz[0])); 00672 MBERRCHK(rval, mb); 00673 00674 } 00675 00676 iBase_EntitySetHandle entityset; //this entityset is for storing the inner nodes on the target surface 00677 vector<iBase_EntityHandle> newNodehandle; 00678 00679 r_err = mk_core()->irel_pair(irelIndx)->getEntSetRelation(targetSurface, 0, entityset); 00680 if (r_err) //there is no entityset associated with targetSurface; this should be an error 00681 { 00682 m_err = mk_core()->imesh_instance()->createEntSet(1, entityset); 00683 IBERRCHK(m_err, "Trouble create the entity set"); 00684 } 00685 00686 bool condition_harmonic = isConcave(); 00687 bool is_exe_harmonic = false; 00688 #ifdef HAVE_ARMADILLO 00689 //target surface mesh projection based on harmonic mapping 00690 if (condition_harmonic){ 00691 is_exe_harmonic = true; 00692 SurfMeshHarmonic(volume, newNodehandle); 00693 } 00694 #endif 00695 00696 if (!is_exe_harmonic){ 00698 // Based on the transformation in the physical space, project the mesh nodes on the source surface to the target surface 00700 Vector3D sPtsCenter(0.0), tPtsCenter(0.0); 00701 std::vector<Vector3D> sBndNodes(sizeBLayer), tBndNodes(sizeBLayer); 00702 int num_pts = 0; 00703 //calculate the barycenters for the source and target boundaries 00704 for (unsigned int i = 0; i < NodeList.size(); i++) 00705 { 00706 if (NodeList[i].onBoundary) 00707 { 00708 sPtsCenter += NodeList[i].xyz; 00709 tPtsCenter += TVertexList[i].xyz; 00710 sBndNodes[num_pts] = NodeList[i].xyz; 00711 tBndNodes[num_pts] = TVertexList[i].xyz; 00712 num_pts++; 00713 } 00714 } 00715 if (sizeBLayer != num_pts) 00716 MBERRCHK(moab::MB_FAILURE, mb); 00717 sPtsCenter = sPtsCenter / num_pts; 00718 tPtsCenter = tPtsCenter / num_pts; 00719 //done with the barycenter calculation 00720 00721 // compute transformation matrix from source boundary to target boundary 00722 Matrix3D transMatrix; 00723 computeTransformationFromSourceToTarget(sBndNodes, sPtsCenter, tBndNodes, tPtsCenter, transMatrix); 00724 00725 //calculate the inner nodes on the target surface 00726 for (unsigned int i = 0; i < NodeList.size(); i++) 00727 { 00728 if (!NodeList[i].onBoundary) 00729 { 00730 Vector3D xyz; 00731 xyz = transMatrix * (NodeList[i].xyz - 2 * sPtsCenter + tPtsCenter)+sPtsCenter; 00732 00733 iGeom::Error g_err = igeom_inst->getEntClosestPt(targetSurface, xyz[0], xyz[1], xyz[2], TVertexList[i].xyz[0], 00734 TVertexList[i].xyz[1], TVertexList[i].xyz[2]); 00735 IBERRCHK(g_err, "Trouble get the closest point on the targets surface."); 00736 00737 //create the node entities 00738 iMesh::Error m_err = mk_core()->imesh_instance()->createVtx(TVertexList[i].xyz[0], TVertexList[i].xyz[1], 00739 TVertexList[i].xyz[2], TVertexList[i].gVertexHandle); 00740 IBERRCHK(m_err, "Trouble create the node entity."); 00741 TVertexList[i].onBoundary = false; 00742 newNodehandle.push_back(TVertexList[i].gVertexHandle); 00743 } 00744 } 00745 } 00746 //add the inner nodes to the entityset 00747 m_err = mk_core()->imesh_instance()->addEntArrToSet(&newNodehandle[0], newNodehandle.size(), entityset); 00748 IBERRCHK(m_err, "Trouble add an array of nodes to the entityset."); 00749 00750 //until now, all the nodes have been generated on the target surface 00751 00752 // we need to decide the orientation with respect to the faces on the source (which are oriented correctly) 00753 // look at the first 2 nodes in the boundary, in target 00754 int sense_out = 1; 00755 std::vector<moab::EntityHandle> bdyTNodes; 00756 std::vector<int> group_sizes; 00757 target_surf-> boundary(0, bdyTNodes, NULL, &group_sizes);// we do not need the senses, yet. Or do we? 00758 // find the first node and second node in the bLayers[] list corresponding to the top 00759 moab::EntityHandle node1 = bdyTNodes[1], node2 = bdyTNodes[2]; 00760 int index = -1; 00761 for (index = 0; index < sizeBLayer; index++) 00762 { 00763 if (bLayers[numLayers * sizeBLayer + index] == node1) 00764 break; 00765 } 00766 if (index == sizeBLayer) 00767 MBERRCHK(moab::MB_FAILURE, mb); 00768 00769 int prevIndex = (index + sizeBLayer -1) % sizeBLayer; 00770 int nextIndex = (index + 1) % sizeBLayer; 00771 if (bLayers[numLayers * sizeBLayer + prevIndex] == node2) 00772 sense_out = -1; 00773 else if (bLayers[numLayers * sizeBLayer + nextIndex] == node2) 00774 sense_out = 1; 00775 else 00776 MBERRCHK(moab::MB_FAILURE, mb);// serious error in the logic, can't decide orientation 00777 00778 00779 //create the quadrilateral elements on the target surface 00780 vector<iBase_EntityHandle> mFaceHandle(FaceList.size()); 00781 vector<iBase_EntityHandle> connect(FaceList.size() * 4); 00782 for (unsigned int i = 0; i < FaceList.size(); i++) 00783 { 00784 if (sense_out < 0) 00785 { 00786 connect[4 * i + 0] = TVertexList[(FaceList[i].getVertex(0))->index].gVertexHandle; 00787 connect[4 * i + 1] = TVertexList[(FaceList[i].getVertex(3))->index].gVertexHandle; 00788 connect[4 * i + 2] = TVertexList[(FaceList[i].getVertex(2))->index].gVertexHandle; 00789 connect[4 * i + 3] = TVertexList[(FaceList[i].getVertex(1))->index].gVertexHandle; 00790 } 00791 else 00792 { 00793 connect[4 * i + 0] = TVertexList[(FaceList[i].getVertex(0))->index].gVertexHandle; 00794 connect[4 * i + 1] = TVertexList[(FaceList[i].getVertex(1))->index].gVertexHandle; 00795 connect[4 * i + 2] = TVertexList[(FaceList[i].getVertex(2))->index].gVertexHandle; 00796 connect[4 * i + 3] = TVertexList[(FaceList[i].getVertex(3))->index].gVertexHandle; 00797 } 00798 00799 } 00800 m_err = mk_core()->imesh_instance()->createEntArr(iMesh_QUADRILATERAL, &connect[0], connect.size(), &mFaceHandle[0]); 00801 IBERRCHK(m_err, "Trouble create the quadrilateral mesh."); 00802 00803 //add the inner face elements to the entityset 00804 m_err = mk_core()->imesh_instance()->addEntArrToSet(&mFaceHandle[0], FaceList.size(), entityset); 00805 IBERRCHK(m_err, "Trouble add an array of quadrilateral entities to the entity set."); 00806 00807 mk_core()->save_mesh("BeforeHex.vtk"); 00808 #ifdef HAVE_MESQUITE 00809 00810 SurfMeshOptimization(); 00811 00812 #endif 00813 00814 return 1; 00815 } 00816 // input: list of nodes on source, boundary center, list of nodes on target, target center 00817 // output: 3x3 matrix A such that 00818 // target= A * ( source - 2*sc + tc) + sc 00819 void OneToOneSwept::computeTransformationFromSourceToTarget(std::vector<Vector3D> & sNodes, Vector3D & sc, 00820 std::vector<Vector3D> & tNodes, Vector3D & tc, Matrix3D & transMatrix) 00821 { 00822 Matrix3D tmpMatrix(0.0); 00823 Matrix3D bMatrix(0.0); 00824 //transform the coordinates 00825 int size = (int)sNodes.size(); 00826 assert(sNodes.size() == tNodes.size()); 00827 for (int i = 0; i < size; i++) 00828 { 00829 //transform the boundary nodes 00830 sNodes[i] = sNodes[i] - 2 * sc + tc; 00831 tNodes[i] = tNodes[i] - sc; 00832 } 00833 00834 //calculate the transformation matrix: transform the nodes on the source surface to the target surface in the physical space 00835 for (int i = 0; i < size; i++) 00836 { 00837 tmpMatrix = tmpMatrix + sNodes[i]*transpose(sNodes[i]); 00838 bMatrix = bMatrix + tNodes[i]*transpose(sNodes[i]); 00839 } 00840 00841 //first determine the affine mapping matrix is singular or not 00842 double detValue = det(tmpMatrix); 00843 assert(pow(detValue, 2)>1.0e-20); 00844 00845 //solve the affine mapping matrix, make use of inverse matrix to get affine mapping matrix 00846 Matrix3D InvMatrix = inverse(tmpMatrix); 00847 transMatrix = bMatrix*InvMatrix; 00848 00849 } 00850 // use similar code to TargetSurfProjection, but do not project on surface... 00851 int OneToOneSwept::ProjectInteriorLayers(std::vector<moab::EntityHandle> & boundLayers, vector<vector<Vertex> > & linkVertexList) 00852 { 00853 iMesh::Error m_err; 00854 Vector3D sPtsCenter(0.), tPtsCenter(0.); 00855 std::vector<Vector3D> PtsCenter(numLayers - 1); 00856 std::vector<Vector3D> sBoundaryNodes(0), tBoundaryNodes(0); 00857 00858 iBase_TagHandle taghandle = 0; 00859 m_err = mk_core()->imesh_instance()->getTagHandle("source", taghandle); 00860 IBERRCHK(m_err, "Trouble getting the source tag handle"); 00861 00862 std::vector<Vector3D> latCoords; 00863 latCoords.resize(boundLayers.size() - 2 * sizeBLayer);// to store the coordinates of the lateral points 00864 // (exclude the source and target, they are already in NodeList and TVertexList)... 00865 00866 moab::ErrorCode rval = mb->get_coords(&(boundLayers[sizeBLayer]), latCoords.size(), &(latCoords[0][0])); //these are the nodes for lateral sides 00867 MBERRCHK(rval, mb); 00868 00869 //calculate the center coordinates 00870 for (int i = 0; i < numLayers - 1; i++) 00871 PtsCenter[i].set(0.); 00872 00873 int numPts = sizeBLayer;// a lot of repetition ; do we really need 00874 sBoundaryNodes.resize(numPts);// another useless copy 00875 tBoundaryNodes.resize(numPts);// another useless copy 00876 for (int j = 0; j < sizeBLayer; j++) 00877 { 00878 iBase_EntityHandle sBNode = (iBase_EntityHandle) boundLayers[j]; 00879 int i; 00880 // this is the index in NodeList... 00881 m_err = mk_core()->imesh_instance()->getIntData(sBNode, taghandle, i); 00882 IBERRCHK(m_err, "Trouble get the int value for nodes in the mesh instance."); 00883 00884 sPtsCenter += NodeList[i].xyz; 00885 tPtsCenter += TVertexList[i].xyz; 00886 00887 sBoundaryNodes[j] = NodeList[i].xyz; 00888 tBoundaryNodes[j] = TVertexList[i].xyz; 00889 00890 for (int lay = 0; lay < numLayers - 1; lay++) 00891 { 00892 // interior layers 00893 int indexNodeOnSide = lay * sizeBLayer + j;// all p3d will be above index i 00894 Vector3D & p3d = latCoords[indexNodeOnSide]; 00895 PtsCenter[lay] += p3d; 00896 linkVertexList[lay][i].xyz = p3d; 00897 linkVertexList[lay][i].gVertexHandle = (iBase_EntityHandle) boundLayers[indexNodeOnSide + sizeBLayer]; 00898 } 00899 } 00900 00901 sPtsCenter = sPtsCenter / numPts; 00902 tPtsCenter = tPtsCenter / numPts; 00903 00904 //calculate the center coordinates for the ith layer 00905 for (int i = 0; i < numLayers - 1; i++) 00906 PtsCenter[i] = PtsCenter[i] / numPts; 00907 00908 std::vector<Vector3D> sNodes(numPts); //boundary nodes on the source surface 00909 std::vector<Vector3D> tNodes(numPts); //boundary nodes on the target surface 00910 std::vector<Vector3D> isBNodes(numPts), itBNodes(numPts); //boundary nodes on the different layers 00911 00912 //loop over different layers 00913 for (int i = 0; i < numLayers - 1; i++) 00914 { 00915 for (int j = 0; j < numPts; j++) 00916 { 00917 sNodes[j] = sBoundaryNodes[j]; 00918 tNodes[j] = tBoundaryNodes[j]; 00919 //coordinates on the different layers 00920 isBNodes[j] = itBNodes[j] = latCoords[i*sizeBLayer + j]; 00921 } 00922 00923 Matrix3D sA, tA; 00924 // from source to layer i 00925 computeTransformationFromSourceToTarget(sNodes, sPtsCenter, isBNodes, PtsCenter[i], sA); 00926 // from target to layer i 00927 computeTransformationFromSourceToTarget(tNodes, tPtsCenter, itBNodes, PtsCenter[i], tA); 00928 //calculate the inner nodes for different layers 00929 double s = (i + 1) / double(numLayers); // interpolation factor 00930 for (unsigned int j = 0; j < NodeList.size(); j++) 00931 { 00932 if (!(NodeList[j].onBoundary)) 00933 { 00934 Vector3D spts, tpts, pts; 00935 iBase_EntityHandle nodeHandle; 00936 spts = sA * (NodeList[j].xyz - 2 * sPtsCenter + PtsCenter[i])+sPtsCenter; 00937 tpts = tA * (TVertexList[j].xyz - 2 * tPtsCenter + PtsCenter[i])+tPtsCenter; 00938 // interpolate the 2 results 00939 pts = (1-s) * spts + s * tpts; 00940 00941 linkVertexList[i][j].xyz = pts; 00942 00943 m_err = mk_core()->imesh_instance()->createVtx(pts[0], pts[1], pts[2], nodeHandle); 00944 IBERRCHK(m_err, "Trouble create the vertex entity."); 00945 linkVertexList[i][j].gVertexHandle = nodeHandle; 00946 m_err = mk_core()->imesh_instance()->addEntToSet(nodeHandle, volumeSet); 00947 IBERRCHK(m_err, "Trouble add the node handle to the entity set."); 00948 } 00949 } 00950 } 00951 return 1; 00952 } 00953 00954 //****************************************************************************// 00955 // function : CreateElements 00956 // Author : Shengyong Cai 00957 // Date : Feb 16, 2011 00958 // Description: create hexahedral elements by connecting 8 nodes which have 00959 // already been created by previous functions 00960 //***************************************************************************// 00961 int OneToOneSwept::CreateElements(vector<vector<Vertex> > &linkVertexList) 00962 { 00963 //create the quadrilateral elements on the different layers 00964 //it is not necessary to create the quadrilateral elements on the different layers. Hex element can be created based on the eight nodes 00965 iMesh::Error m_err; 00966 vector<iBase_EntityHandle> mVolumeHandle(FaceList.size()); 00967 // first decide orientation, based on the FaceList orientation, and first node above 00968 // take one face on source, and first node above in layer 1 00969 Vertex * v1 = FaceList[0].connect[0]; 00970 Vertex * v2 = FaceList[0].connect[1]; 00971 Vertex * v3 = FaceList[0].connect[2]; 00972 // vertex 4 is on layer 1, above vertex 1 00973 int ind1 = v1->index; 00974 Vertex & v4 = linkVertexList[0][ind1];// this is the node above vertex 1 in layer 1 00975 00976 Vector3D normal1 = (v2->xyz-v1->xyz)*(v3->xyz-v1->xyz); 00977 double vol6= normal1 % (v4.xyz-v1->xyz); 00978 std::cout << "Orientation is "; 00979 int skip = 0; 00980 if (vol6 < 0) 00981 { 00982 std::cout <<"negative\n"; 00983 skip=4; 00984 } 00985 else 00986 std::cout <<"positive\n"; 00987 00988 vector<iBase_EntityHandle> connect(8); 00989 for (int m = 0; m < numLayers; m++) 00990 { 00991 if (m == 0) // first layer, above source 00992 { 00993 for (unsigned int i = 0; i < FaceList.size(); i++) 00994 { 00995 for (int k=0; k<4; k++) 00996 { 00997 connect[k+skip]=NodeList[(FaceList[i].connect[k])->index].gVertexHandle; 00998 connect[(k+skip+4)%8] = linkVertexList[m][(FaceList[i].connect[k])->index].gVertexHandle; 00999 } 01000 m_err = mk_core()->imesh_instance()->createEnt(iMesh_HEXAHEDRON, &connect[0], 8, mVolumeHandle[i]); 01001 IBERRCHK(m_err, "Trouble create the hexahedral element."); 01002 } 01003 m_err = mk_core()->imesh_instance()->addEntArrToSet(&mVolumeHandle[0], FaceList.size(), volumeSet); 01004 IBERRCHK(m_err, "Trouble add an array of hexahedral elements to the entity set."); 01005 } 01006 else if (m == (numLayers - 1)) // the top layer, below target 01007 { 01008 for (unsigned int i = 0; i < FaceList.size(); i++) 01009 { 01010 for (int k=0; k<4; k++) 01011 { 01012 connect[k+skip]=linkVertexList[m-1][(FaceList[i].connect[k])->index].gVertexHandle; 01013 connect[(k+skip+4)%8] = TVertexList[(FaceList[i].connect[k])->index].gVertexHandle; 01014 } 01015 m_err = mk_core()->imesh_instance()->createEnt(iMesh_HEXAHEDRON, &connect[0], 8, mVolumeHandle[i]); 01016 IBERRCHK(m_err, "Trouble create the hexahedral elements."); 01017 } 01018 m_err = mk_core()->imesh_instance()->addEntArrToSet(&mVolumeHandle[0], FaceList.size(), volumeSet); 01019 IBERRCHK(m_err, "Trouble add an array of hexahedral elements to the entity set."); 01020 } 01021 else // intermediate layers, m is between 0 and num_layers-2 01022 { 01023 for (unsigned int i = 0; i < FaceList.size(); i++) 01024 { 01025 for (int k=0; k<4; k++) 01026 { 01027 connect[k+skip]=linkVertexList[m-1][(FaceList[i].connect[k])->index].gVertexHandle; 01028 connect[(k+skip+4)%8] = linkVertexList[m][(FaceList[i].connect[k])->index].gVertexHandle; 01029 } 01030 m_err = mk_core()->imesh_instance()->createEnt(iMesh_HEXAHEDRON, &connect[0], 8, mVolumeHandle[i]); 01031 IBERRCHK(m_err, "Trouble create the hexahedral elements."); 01032 } 01033 m_err = mk_core()->imesh_instance()->addEntArrToSet(&mVolumeHandle[0], FaceList.size(), volumeSet); 01034 IBERRCHK(m_err, "Trouble add an array of hexahedral elements to the entity set."); 01035 } 01036 } 01037 return 1; 01038 } 01039 01040 #ifdef HAVE_MESQUITE 01041 //target surface mesh smoothing by Mesquite 01042 void OneToOneSwept::SurfMeshOptimization() 01043 { 01044 MEntSelection::iterator mit = mentSelection.begin(); 01045 ModelEnt *me = mit -> first; 01046 int irelPairIndex = me->iRelPairIndex(); 01047 iGeom * igeom_inst = mk_core()->igeom_instance(me->iGeomIndex()); 01048 //create a tag to attach the coordinates to nodes 01049 iBase_TagHandle mapped_tag = 0; 01050 iMesh::Error m_err = mk_core()->imesh_instance()->getTagHandle("MsqAltCoords", mapped_tag); 01051 if (m_err) 01052 { 01053 m_err = mk_core()->imesh_instance()->createTag("MsqAltCoords", 3, iBase_DOUBLE, mapped_tag); 01054 IBERRCHK(m_err, "Trouble create a tag."); 01055 01056 } 01057 //attach the coordinates to the nodes 01058 double tag_data[3*TVertexList.size()]; 01059 std::vector<iBase_EntityHandle> vertexHandle(TVertexList.size()); 01060 for (unsigned int i = 0; i < NodeList.size();i++) 01061 { 01062 tag_data[3*i] = NodeList[i].xyz[0]; 01063 tag_data[3*i+1] = NodeList[i].xyz[1]; 01064 tag_data[3*i+2] = NodeList[i].xyz[2]; 01065 vertexHandle[i] = TVertexList[i].gVertexHandle; 01066 } 01067 m_err = mk_core()->imesh_instance()->setDblArrData(&vertexHandle[0], NodeList.size(), mapped_tag, &tag_data[0]); 01068 IBERRCHK(m_err, "Trouble set an array of int data to nodes."); 01069 01070 //get the mesh entityset for target surface 01071 iBase_EntitySetHandle surfSets; 01072 iRel::Error r_err = mk_core()->irel_pair(irelPairIndex)->getEntSetRelation(targetSurface, 0, surfSets); 01073 IBERRCHK(r_err, "Trouble get the mesh entity set for the target surface."); 01074 //call the MeshImprove class to smooth the target surface mesh by using Mesquite 01075 MeshImprove meshopt(mk_core(), false, true, true, true, igeom_inst); 01076 meshopt.SurfMeshImprove(targetSurface, surfSets, iBase_FACE); 01077 } 01078 #endif 01079 } 01080