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