MeshKit  1.0
OneToOneSwept.cpp
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines