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