MeshKit  1.0
example_fbgeom.cpp
Go to the documentation of this file.
00001 
00014 #include "meshkit/MKCore.hpp"
00015 #include "meshkit/FBiGeom.hpp"
00016 #include "meshkit/SizingFunction.hpp"
00017 #include "meshkit/ModelEnt.hpp"
00018 
00019 using namespace MeshKit;
00020 
00021 #include "stdlib.h"
00022 #include <iostream>
00023 #include <sstream>
00024 #include <fstream>
00025 #include <math.h>
00026 #include <assert.h>
00027 #include <string.h>
00028 #include <stdio.h>
00029 #include <time.h>
00030 #include "moab/CartVect.hpp"
00031 #include "moab/Interface.hpp"
00032 #include "moab/Skinner.hpp"
00033 //#include "MBiMesh.hpp"
00034 #include <vector>
00035 
00036 MKCore * mk;
00037 bool save_mesh = true;
00038 //std::string file_name;
00039 std::string bottom_mesh; //="quadsOnBed.vtk";
00040 std::string top_filename; // "cropSE.h5m", to be smaller
00041 std::string out_mesh; // deleted by default
00042 int nbLayers;
00043 std::vector<double> grades;
00044 
00045 bool create_a_set_with_tag_value(iMesh_Instance meshIface,
00046                                  iBase_EntityHandle * ents, int num_ents, char * tag_name,
00047                                  int size_name_tag, int value, iBase_EntitySetHandle & mesh_set);
00048 
00049 int main(int argc, char *argv[])
00050 {
00051     // Check command line arg
00052 
00053     nbLayers = 5; // a good number of layers
00054     if (argc < 4) {
00055         std::cout << "Usage: " << argv[0]
00056                   << " <bottom_mesh> <top_filename> <out_mesh> [-n layers]  <ratios> "
00057 
00058                   << std::endl;
00059         bottom_mesh = (std::string) MESH_DIR + "/" + "quadsOnBed.vtk";
00060         top_filename = std::string (MESH_DIR) + "/" + "cropSE.h5m"; //, to be smaller
00061         out_mesh = "hexas.h5m"; // deleted by default
00062         std::cout << " Using : " << bottom_mesh << " " << top_filename << " "
00063                   << out_mesh << " -n 5" << "\n";
00064 
00065         std::cout << " the default output file is deleted \n";
00066     } else {
00067         bottom_mesh = argv[1];
00068         top_filename = argv[2];
00069         out_mesh = argv[3];
00070         save_mesh = true;
00071         int argno = 4;
00072         while (argno < argc) {
00073             if (!strcmp(argv[argno], "-n")) {
00074                 argno++;
00075                 nbLayers = atoi(argv[argno]);
00076                 argno++;
00077             } else {
00078                 //everything that is not interpreted must be a ratio
00079                 grades.push_back(atof(argv[argno]));
00080                 argno++;
00081             }
00082         }
00083     }
00084     int nbGrades = (int) grades.size();
00085     if (nbGrades > 0) {
00086         // add all grades, then divide by total to get the new accumulated ratios
00087         double total = 0.;
00088         int i;
00089         for (i = 0; i < nbGrades; i++) {
00090             total += grades[i];
00091             //grades[i] = total; // accumulated so far
00092         }
00093         for (i = 0; i < nbGrades; i++) {
00094             grades[i] /= total;
00095         }
00096         nbLayers = nbGrades;
00097     } else {
00098         nbGrades = nbLayers;
00099         for (int i = 0; i < nbLayers; i++)
00100             grades.push_back(1. / nbLayers);
00101     }
00102 
00103     clock_t start_time = clock();
00104     mk = new MKCore;
00105     iBase_ErrorType err;
00106     FBiGeom * fbiGeom = new FBiGeom(); // true for smooth, false for linear
00107     mk->add_igeom_instance(fbiGeom);
00108     // read initial mesh (quad mesh surface bottom)
00109     char opts[]="SMOOTH;";
00110     err = fbiGeom->load(top_filename.c_str(), opts);// loading the top surface , as geometry
00111     IBERRCHK(err, "Failure to load smooth geometry.");
00112     clock_t load_time1 = clock();// load fb file (expensive),
00113     // load here the quads from bottom
00114     iBase_EntitySetHandle setForQuads;
00115 
00116     err = mk->imesh_instance()->createEntSet(false, setForQuads);
00117     IBERRCHK(err, "Failure to create set for loading the quads");
00118     // load the quads in this set
00119     err = mk->imesh_instance()->load(setForQuads, bottom_mesh.c_str());
00120     IBERRCHK(err, "Failure to load quad file.");
00121 
00122     clock_t load_time2 = clock();
00123     // and quads file (cheap)
00124     //
00125     double direction[3] = { 0., 0., 1. }; // normalized
00126 
00127     // can we get to iMesh from iMesh.hpp?
00128     iMesh_Instance mesh1 = mk->imesh_instance()->instance();
00129 
00130     // first get all nodes from bottom mesh
00131     // get the quads and the vertices in one shot
00132     iBase_EntityHandle *quads = NULL;
00133     int quads_alloc = 0;
00134     iBase_EntityHandle *vert_adj = NULL;
00135     int vert_adj_alloc = 0, vert_adj_size, numQuads;
00136     int * offsets = NULL, offsets_alloc = 0, indices_size;
00137     int * indices = NULL, indices_alloc = 0, offsets_size;
00138     int ierr=0;
00139     iMesh_getAdjEntIndices(mesh1, setForQuads, iBase_FACE, iMesh_QUADRILATERAL,
00140                            iBase_VERTEX, &quads, &quads_alloc, &numQuads, &vert_adj,
00141                            &vert_adj_alloc, &vert_adj_size, &indices, &indices_alloc, &indices_size,
00142                            &offsets, &offsets_alloc, &offsets_size, &ierr);
00143     // create one set with desired tag (1)
00144     iBase_EntitySetHandle mesh_set1;
00145     char nsname[]= "NEUMANN_SET";
00146     create_a_set_with_tag_value(mesh1, quads, numQuads, nsname , 11, 1, mesh_set1);
00147     /* get the coordinates in one array */
00148 
00149     // delete now the original set?
00150     int vert_coords_alloc = 0;
00151     double * xyz = 0; // not allocated
00152     int vertex_coord_size = 0;
00153 
00154     iMesh_getVtxArrCoords(mesh1, vert_adj, vert_adj_size, iBase_INTERLEAVED,
00155                           &xyz, &vert_coords_alloc, &vertex_coord_size, &ierr);
00156     //IBERRCHK(err, "failed to get vertex coordinates of entities.");
00157 
00158     // then, go to shoot rays to decide the sweeping thickness
00159     int & numNodes = vert_adj_size; // to not change too much
00160     double * dArr = new double[numNodes];
00161     // then, go to shoot rays
00162     int j = 0;
00163     int numRaysIntersected = 0;
00164     // double factorFloating = (1.-937./1026.);
00165     for (j = 0; j < numNodes; j++) {
00166         //dArr[j] = xyz[j * 3 + 2];
00167         dArr[j] = 0; // no intersection situation, marked by 0
00168         // for a point, see if it is inside the polygon, with winding number
00169 
00170         iBase_StorageOrder order=iBase_INTERLEAVED;
00171         std::vector<iGeom::EntityHandle> entities_out;
00172         std::vector<double> points_out;
00173         std::vector<double> params_out;
00174 
00175         double pos[3] = { xyz[3 * j], xyz[3 * j + 1], xyz[3 * j + 2] };
00176         err= fbiGeom->getPntRayIntsct(pos[0], pos[1], pos[2], direction[0],
00177                 direction[1], direction[2], order, entities_out, points_out,
00178                 params_out);
00179         // get the first coordinate
00180         if (err != 0 || entities_out.size() < 1)
00181             continue;// we should bail out
00182         //double zTop = points_out[2]; // the z of the top, intersection computation
00183 
00184         // consider only the first intersection point
00185         dArr[j] = params_out[0]; // the first intersection only
00186         numRaysIntersected++;
00187     }
00188     // to xyz coordinates, add some thickness, and (thick/layers), to create
00189     // a next layer, and the hexas
00190     iBase_EntityHandle * layer1 = vert_adj;
00191     iBase_EntitySetHandle setForHexas;
00192 
00193     err = mk->imesh_instance()->createEntSet(false, setForHexas);
00194     IBERRCHK(err, "Failure to create set for loading the hexas");
00195     // create in a loop layer 2 , and vertices in layer 2, and hexas
00196     // between layer 1 and layer2
00197     iBase_EntitySetHandle mesh_set2;
00198     for (int i = 0; i < nbLayers; i++) {
00199         for (int j = 0; j < numNodes; j++) {
00200             for (int k = 0; k < 3; k++)
00201                 xyz[3 * j + k] += direction[k] * dArr[j] * grades[i];
00202         }
00203         // now create new vertices at this position, layer 2
00204         iBase_EntityHandle * newVerts = NULL; // no vertices yet
00205         int size1, size2;
00206         iMesh_createVtxArr(mesh1,
00207                            /*in*/numNodes,
00208                            /*in*/iBase_INTERLEAVED,
00209                            /*in*/xyz,
00210                            /*in*/numNodes * 3,
00211                            /*inout*/&newVerts,
00212                            /*inout*/&size1,
00213                            /*inout*/&size2,
00214                            /*out*/&ierr);
00215 
00216         int numHexa = numQuads;
00217         long int * adjacency = new long int[8 * numHexa];
00218         iBase_EntityHandle * conn = (iBase_EntityHandle *) adjacency;
00219         for (int L = 0; L < numHexa; L++) {
00220 
00221             for (int k = 0; k < 4; k++) {
00222                 conn[8 * L + k] = layer1[indices[offsets[L] + k]];
00223                 conn[8 * L + k + 4] = newVerts[indices[offsets[L] + k]];
00224             }
00225         }
00226 
00227         int n = numHexa;
00228         int junk1 = n, junk2 = n, junk3 = n, junk4 = n;
00229         int * stat = new int[numHexa];
00230         int* ptr2 = stat;
00231         //int ierr;
00232         iBase_EntityHandle * new_entity_handles = NULL;
00233         iMesh_createEntArr(mesh1, iMesh_HEXAHEDRON, conn, 8 * n,
00234                            &new_entity_handles, &junk1, &junk2, &ptr2, &junk3, &junk4, &ierr);
00235         // add hexas to the setForHexas set
00236         iMesh_addEntArrToSet(mesh1,
00237                              new_entity_handles, numHexa,
00238                              setForHexas,
00239                              &ierr );
00240         // end copy
00241         // at the end, layer 1 becomes layer 2
00242         layer1 = newVerts;
00243         // if last layer, create top quads too
00244         if (i == nbLayers - 1) {
00245             // last layer, create top quads too
00246             for (int L = 0; L < numQuads; L++) {
00247 
00248                 for (int k = 0; k < 4; k++) {
00249                     //conn[8*L+k]   = layer1[ indices[offsets[L]+k ] ];
00250                     conn[4 * L + k] = newVerts[indices[offsets[L] + k]];
00251                 }
00252             }
00253 
00254             iMesh_createEntArr(mesh1, iMesh_QUADRILATERAL, conn, 4 * n,
00255                                &new_entity_handles, &junk1, &junk2, &ptr2, &junk3, &junk4, &ierr);
00256             // create another set
00257             create_a_set_with_tag_value(mesh1, new_entity_handles, numQuads,
00258                                         nsname, 11, 2, mesh_set2); // top
00259         }
00260     }
00261 
00262     clock_t compute_time = clock();
00263 
00264     // now get MOAB , to use skin on the hexas, to get the rest of quads
00265     moab::Interface * mb = mk->moab_instance();
00266     // get all hexas, and get the skin
00267     moab::Range hexas, iniQuads;
00268     mb->get_entities_by_type(0, moab::MBHEX, hexas);
00269     mb->get_entities_by_type(0, moab::MBQUAD, iniQuads);
00270     moab::Skinner skinner(mb);
00271     moab::Range allQuads;
00272     skinner.find_skin(0, hexas, 2, allQuads);
00273 
00274     moab::Range latQuads = subtract(allQuads, iniQuads);
00275 
00276     // add this one too to the hexa set
00277     moab::EntityHandle latSet;
00278     mb->create_meshset(moab::MESHSET_SET, latSet);
00279     mb->add_entities(latSet, latQuads);
00280 
00281     // add range and tag
00282     moab::Tag ntag;
00283     mb->tag_get_handle("NEUMANN_SET", 1,
00284                        moab::MB_TYPE_INTEGER,  ntag);
00285     int val = 3;
00286     mb->tag_set_data(ntag, &latSet, 1, (void*) &val);
00287 
00288     // to the hexa set, add chlidren sets, mesh_set1, etc
00289     mb->add_parent_child( (moab::EntityHandle)setForHexas,
00290                           (moab::EntityHandle) mesh_set1);
00291     mb->add_parent_child( (moab::EntityHandle)setForHexas,
00292                           (moab::EntityHandle) mesh_set2);
00293     mb->add_parent_child( (moab::EntityHandle)setForHexas,
00294                           latSet);
00295 
00296     delete fbiGeom;
00297     assert(iBase_SUCCESS == err);
00298     iMesh_save(mesh1, setForHexas, out_mesh.c_str(), 0, &ierr, out_mesh.length(), 0);
00299     if (iBase_SUCCESS != err) {
00300         std::cerr << "ERROR saving mesh to " << out_mesh << std::endl;
00301         return 0;
00302     }
00303     clock_t out_time = clock();
00304     std::cout << "Total time is " << (double) (out_time - start_time)
00305                  / CLOCKS_PER_SEC << " s\n  load top : " << (double) (load_time1
00306                                                                       - start_time) / CLOCKS_PER_SEC << " s\n  load bottom : "
00307               << (double) (load_time2 - load_time1) / CLOCKS_PER_SEC
00308               << " s\n  compute time : " << (double) (compute_time - load_time2)
00309                  / CLOCKS_PER_SEC << " s\n  write time : " << (double) (out_time
00310                                                                         - compute_time) / CLOCKS_PER_SEC << std::endl;
00311     std::cout << "num rays intersected : " << numRaysIntersected << "\n";
00312 
00313     free(xyz);
00314     delete[] dArr;
00315 
00316     if (!save_mesh)
00317     {
00318         remove(out_mesh.c_str());
00319     }
00320     return 0;
00321 }
00322 
00323 // Helper function: usually a NEUMANN set
00324 bool create_a_set_with_tag_value(iMesh_Instance meshIface,
00325                                  iBase_EntityHandle * ents, int num_ents, char * tag_name,
00326                                  int size_name_tag, int value, iBase_EntitySetHandle & mesh_set)
00327 {
00328     // will first create a new entity set in the mesh; will add the
00329     // entities;
00330     // will associate the tag with the given value (create the tag if not existent yet)
00331     // NEUMANN_SET
00332     int result = iBase_SUCCESS;
00333 
00334     bool isList = false;
00335     //iBase_EntitySetHandle mesh_set;
00336 
00337     iMesh_createEntSet(meshIface, isList, &mesh_set, &result);
00338     if (iBase_SUCCESS != result)
00339         return false;
00340 
00341     iBase_TagHandle tag_handle;
00342     iMesh_getTagHandle(meshIface, tag_name, &tag_handle, &result, size_name_tag);
00343     assert (0 != tag_handle);
00344     iMesh_setEntSetIntData(meshIface, mesh_set, tag_handle, value, &result);
00345     if (iBase_SUCCESS != result)
00346         return false;
00347     // add the entities to the set
00348     //
00349     iMesh_addEntArrToSet(meshIface, ents, num_ents, mesh_set, &result);
00350     if (iBase_SUCCESS != result)
00351         return false;
00352 
00353     return true;
00354 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines