moab
ReadSms.cpp
Go to the documentation of this file.
00001 
00022 #include "ReadSms.hpp"
00023 #include "FileTokenizer.hpp" // for file tokenizer
00024 #include "Internals.hpp"
00025 #include "moab/Interface.hpp"
00026 #include "moab/ReadUtilIface.hpp"
00027 #include "moab/Range.hpp"
00028 #include "MBTagConventions.hpp"
00029 #include "MBParallelConventions.h"
00030 #include "moab/CN.hpp"
00031 
00032 #include <errno.h>
00033 #include <string.h>
00034 #include <map>
00035 #include <set>
00036 #include <iostream>
00037 
00038 #define CHECK(a) if (MB_SUCCESS != result) {      \
00039       std::cerr << a << std::endl;                \
00040       return result;                              \
00041     }
00042 
00043 #define CHECKN(a) if (n != (a)) return MB_FILE_WRITE_ERROR
00044 
00045 namespace moab {
00046 
00047 ReaderIface* ReadSms::factory( Interface* iface )
00048   { return new ReadSms(iface); }
00049 
00050 ReadSms::ReadSms(Interface* impl)
00051     : mdbImpl(impl)
00052 {
00053   mdbImpl->query_interface(readMeshIface);
00054 }
00055 
00056 ReadSms::~ReadSms()
00057 {
00058   if (readMeshIface) {
00059     mdbImpl->release_interface(readMeshIface);
00060     readMeshIface = 0;
00061   }
00062 }
00063 
00064 ErrorCode ReadSms::read_tag_values(const char* /* file_name */,
00065                                    const char* /* tag_name */,
00066                                    const FileOptions& /* opts */,
00067                                    std::vector<int>& /* tag_values_out */,
00068                                    const SubsetList* /* subset_list */ )
00069 {
00070   return MB_NOT_IMPLEMENTED;
00071 }
00072 
00073 
00074 ErrorCode ReadSms::load_file( const char* filename, 
00075                               const EntityHandle* ,
00076                               const FileOptions& ,
00077                               const ReaderIface::SubsetList* subset_list,
00078                               const Tag* file_id_tag )
00079 {
00080   if (subset_list) {
00081     readMeshIface->report_error( "Reading subset of files not supported for Sms." );
00082     return MB_UNSUPPORTED_OPERATION;
00083   }
00084 
00085   setId = 1;
00086     
00087     // Open file
00088   FILE* file_ptr = fopen( filename, "r" );
00089   if (!file_ptr)
00090   {
00091     readMeshIface->report_error( "%s: %s\n", filename, strerror(errno) );
00092     return MB_FILE_DOES_NOT_EXIST;
00093   }
00094 
00095   const ErrorCode result = load_file_impl( file_ptr, file_id_tag );
00096   fclose( file_ptr );
00097 
00098   return result;
00099 }
00100 
00101 ErrorCode ReadSms::load_file_impl( FILE* file_ptr, const Tag* file_id_tag )
00102 {
00103   bool warned = false;
00104   double dum_params[] = {0.0, 0.0, 0.0};
00105   
00106   int zero = 0;
00107   ErrorCode result = mdbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER,
00108                                               globalId, MB_TAG_DENSE|MB_TAG_CREAT, &zero);
00109   CHECK("Failed to create gid tag.");
00110     
00111   result = mdbImpl->tag_get_handle("PARAMETER_COORDS", 3, MB_TYPE_DOUBLE,
00112                                    paramCoords, MB_TAG_DENSE|MB_TAG_CREAT );
00113   CHECK("Failed to create param coords tag.");
00114     
00115   int negone = -1;
00116   result = mdbImpl->tag_get_handle(GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER,
00117                                    geomDimension, MB_TAG_SPARSE|MB_TAG_CREAT, &negone);
00118   CHECK("Failed to create geom dim tag.");
00119 
00120   int n;
00121   char line[256];
00122   int file_type;
00123   n = fscanf(file_ptr, "%s %d", line, &file_type);
00124   CHECKN(2);
00125 
00126   if (3 == file_type) {
00127     result = read_parallel_info(file_ptr);
00128     CHECK("Failed to read parallel info.");
00129   }
00130 
00131   int nregions, nfaces, nedges, nvertices, npoints;
00132   n = fscanf(file_ptr, "%d %d %d %d %d", &nregions, &nfaces, &nedges,
00133          &nvertices, &npoints);
00134   CHECKN(5);
00135   if (nregions < 0 || nfaces < 0 || nedges < 0 || nvertices < 0 || npoints < 0)
00136     return MB_FILE_WRITE_ERROR;
00137 
00138     // create the vertices
00139   std::vector<double*> coord_arrays;
00140   EntityHandle vstart = 0;
00141   result = readMeshIface->get_node_coords( 3, nvertices, MB_START_ID, 
00142                                            vstart, coord_arrays );
00143   CHECK("Failed to get node arrays.");
00144   if (MB_SUCCESS != result)
00145     return result;
00146     
00147   if (file_id_tag) {
00148     result = add_entities( vstart, nvertices, file_id_tag );
00149     if (MB_SUCCESS != result)
00150       return result;
00151   }
00152   
00153   EntityHandle this_gent, new_handle;
00154   std::vector<EntityHandle> gentities[4];
00155   int gent_id, dum_int;
00156   int gent_type, num_connections;
00157   
00158   for(int i = 0; i < nvertices; i++)
00159   {
00160     n = fscanf(file_ptr, "%d", &gent_id); 
00161     CHECKN(1);
00162     if (!gent_id) continue;
00163 
00164     n = fscanf(file_ptr,"%d %d %lf %lf %lf", &gent_type, &num_connections,
00165            coord_arrays[0]+i, coord_arrays[1]+i, coord_arrays[2]+i);
00166     CHECKN(5);
00167     
00168     result = get_set(gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );
00169     if (MB_SUCCESS != result)
00170       return result;
00171 
00172     new_handle = vstart + i;
00173     result = mdbImpl->add_entities(this_gent, &new_handle, 1);
00174     CHECK("Adding vertex to geom set failed.");
00175     if (MB_SUCCESS != result) return result;
00176     
00177     switch(gent_type)
00178     {
00179       case 1:
00180           n = fscanf(file_ptr, "%le", dum_params);
00181           CHECKN(1);
00182           result = mdbImpl->tag_set_data(paramCoords, &new_handle, 1, dum_params);
00183           CHECK("Failed to set param coords tag for vertex.");
00184           if (MB_SUCCESS != result) return result;
00185           break;
00186       case 2:
00187           n = fscanf(file_ptr, "%le %le %d", dum_params, dum_params+1, &dum_int);
00188           CHECKN(3);
00189           dum_params[2] = dum_int;
00190           result = mdbImpl->tag_set_data(paramCoords, &new_handle, 1, dum_params);
00191           CHECK("Failed to set param coords tag for vertex.");
00192           if (MB_SUCCESS != result) return result;
00193           break;
00194           
00195       default: break;
00196     }
00197   } // end of reading vertices
00198 
00199 // *******************************
00200 //  Read Edges
00201 // *******************************
00202 
00203   int vert1, vert2, num_pts;
00204   std::vector<EntityHandle> everts(2);
00205   EntityHandle estart, *connect;
00206   result = readMeshIface->get_element_connect(nedges, 2, MBEDGE, 1, estart, connect);
00207   CHECK("Failed to create array of edges.");
00208   if (MB_SUCCESS != result) return result;
00209 
00210   if (file_id_tag) {
00211     result = add_entities( estart, nedges, file_id_tag );
00212     if (MB_SUCCESS != result)
00213       return result;
00214   }
00215   
00216   for(int i = 0; i < nedges; i++)
00217   {
00218     n = fscanf(file_ptr,"%d",&gent_id);
00219     CHECKN(1);
00220     if (!gent_id) continue;
00221 
00222     n = fscanf(file_ptr, "%d %d %d %d %d", &gent_type, &vert1, &vert2, 
00223            &num_connections, &num_pts);
00224     CHECKN(5);
00225     if (vert1 < 1 || vert1 > nvertices)
00226       return MB_FILE_WRITE_ERROR;
00227     if (vert2 < 1 || vert2 > nvertices)
00228       return MB_FILE_WRITE_ERROR;
00229     
00230     connect[0] = vstart + vert1 - 1;
00231     connect[1] = vstart + vert2 - 1;
00232     if (num_pts > 1 && !warned) {
00233       std::cout << "Warning: num_points > 1 not supported; choosing last one." << std::endl;
00234       warned = true;
00235     }
00236 
00237     result = get_set(gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag);
00238     CHECK("Problem getting geom set for edge.");
00239     if (MB_SUCCESS != result)
00240       return result;
00241 
00242     new_handle = estart + i;
00243     result = mdbImpl->add_entities(this_gent, &new_handle, 1);
00244     CHECK("Failed to add edge to geom set.");
00245     if (MB_SUCCESS != result) return result;
00246 
00247     connect += 2;
00248 
00249     for(int j = 0; j < num_pts; j++) {
00250       switch(gent_type) {
00251         case 1: 
00252             n = fscanf(file_ptr, "%le", dum_params);
00253             CHECKN(1);
00254             result = mdbImpl->tag_set_data(paramCoords, &new_handle, 1, dum_params);
00255             CHECK("Failed to set param coords tag for edge.");
00256             if (MB_SUCCESS != result) return result;
00257             break;
00258         case 2: 
00259             n = fscanf(file_ptr, "%le %le %d", dum_params, dum_params+1, &dum_int);
00260             CHECKN(3);
00261             dum_params[2] = dum_int;
00262             result = mdbImpl->tag_set_data(paramCoords, &new_handle, 1, dum_params);
00263             CHECK("Failed to set param coords tag for edge.");
00264             if (MB_SUCCESS != result) return result;
00265             break;
00266         default: 
00267             break;
00268       }
00269     }
00270 
00271   }  // end of reading edges
00272 
00273 // *******************************
00274 //  Read Faces
00275 // *******************************
00276   std::vector<EntityHandle> bound_ents, bound_verts, new_faces;
00277   int bound_id;
00278   Range shverts;
00279   new_faces.resize(nfaces);
00280   int num_bounding;
00281     
00282   for(int i = 0; i < nfaces; i++)
00283   {
00284     n = fscanf(file_ptr, "%d", &gent_id);
00285     CHECKN(1);
00286     if(!gent_id) continue;
00287 
00288     n = fscanf(file_ptr,"%d %d", &gent_type, &num_bounding);
00289     CHECKN(2);
00290 
00291     result = get_set(gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag);
00292     CHECK("Problem getting geom set for face.");
00293     if (MB_SUCCESS != result)
00294       return result;
00295 
00296     bound_ents.resize(num_bounding+1);
00297     bound_verts.resize(num_bounding);
00298     for(int j = 0; j < num_bounding; j++) {
00299       n = fscanf(file_ptr, "%d ", &bound_id);
00300       CHECKN(1);
00301       if (0 > bound_id) bound_id = abs(bound_id);
00302       assert(0 < bound_id && bound_id <= nedges);
00303       if (bound_id < 1 || bound_id > nedges)
00304         return MB_FILE_WRITE_ERROR;
00305       bound_ents[j] = estart + abs(bound_id) - 1;
00306     }
00307 
00308       // convert edge-based model to vertex-based one
00309     for (int j = 0; j < num_bounding; j++) {
00310       if (j == num_bounding-1) bound_ents[j+1] = bound_ents[0];
00311       result = mdbImpl->get_adjacencies(&bound_ents[j], 2, 0, false, shverts);
00312       CHECK("Failed to get vertices bounding edge.");
00313       if (MB_SUCCESS != result) return result;
00314       assert(shverts.size() == 1);
00315       bound_verts[j] = *shverts.begin();
00316       shverts.clear();
00317     }
00318 
00319     result = mdbImpl->create_element((EntityType)(MBTRI+num_bounding-3),
00320                                      &bound_verts[0], bound_verts.size(), 
00321                                      new_faces[i]);
00322     CHECK("Failed to create edge.");
00323     if (MB_SUCCESS != result) return result;
00324 
00325     result = mdbImpl->add_entities(this_gent, &new_faces[i], 1);
00326     CHECK("Failed to add edge to geom set.");
00327     if (MB_SUCCESS != result) return result;
00328 
00329     int num_read = fscanf(file_ptr, "%d", &num_pts);
00330     if(!num_pts || !num_read) continue;
00331 
00332     for(int j = 0; j < num_pts; j++) {
00333       switch(gent_type) {
00334         case 1: 
00335             n = fscanf(file_ptr, "%le", dum_params); CHECKN(1);
00336             result = mdbImpl->tag_set_data(paramCoords, &new_faces[i], 1, dum_params);
00337             CHECK("Failed to set param coords tag for face.");
00338             if (MB_SUCCESS != result) return result;
00339             break;
00340         case 2: 
00341             n = fscanf(file_ptr, "%le %le %d", dum_params, dum_params+1, &dum_int);
00342             CHECKN(3);
00343             dum_params[2] = dum_int;
00344             result = mdbImpl->tag_set_data(paramCoords, &new_faces[i], 1, dum_params);
00345             CHECK("Failed to set param coords tag for face.");
00346             if (MB_SUCCESS != result) return result;
00347             break;
00348         default: 
00349             break;
00350       }
00351     }
00352 
00353   } // end of reading faces
00354   
00355   if (file_id_tag) {
00356     result = readMeshIface->assign_ids( *file_id_tag, &new_faces[0], new_faces.size(), 1 );
00357     if (MB_SUCCESS != result)
00358       return result;
00359   }
00360 
00361 
00362 // *******************************
00363 //  Read Regions
00364 // *******************************
00365   int sense[MAX_SUB_ENTITIES];
00366   bound_verts.resize(MAX_SUB_ENTITIES);
00367 
00368   std::vector<EntityHandle> regions;
00369   if (file_id_tag)
00370     regions.resize( nregions );
00371   for(int i = 0; i < nregions; i++)
00372   {
00373     n = fscanf(file_ptr, "%d", &gent_id); CHECKN(1);
00374     if (!gent_id) continue;
00375     result = get_set(gentities, 3, gent_id, geomDimension, this_gent, file_id_tag);
00376     CHECK("Couldn't get geom set for region.");
00377     if (MB_SUCCESS != result)
00378       return result;
00379     n = fscanf(file_ptr, "%d", &num_bounding); CHECKN(1);
00380     bound_ents.resize(num_bounding);
00381     for(int j = 0; j < num_bounding; j++) {
00382       n = fscanf(file_ptr, "%d ", &bound_id); CHECKN(1);
00383       assert(abs(bound_id) < (int)new_faces.size()+1 && bound_id);
00384       if (!bound_id || abs(bound_id) > nfaces)
00385         return MB_FILE_WRITE_ERROR;
00386       sense[j] = (bound_id < 0) ? -1 : 1;
00387       bound_ents[j] = new_faces[abs(bound_id)-1];
00388     }
00389 
00390     EntityType etype;
00391     result = readMeshIface->get_ordered_vertices(&bound_ents[0], sense, 
00392                                                  num_bounding,
00393                                                  3, &bound_verts[0], etype);
00394     CHECK("Failed in get_ordered_vertices.");
00395         
00396       // make the element
00397     result = mdbImpl->create_element(etype, &bound_verts[0], 
00398                                      CN::VerticesPerEntity(etype), new_handle);
00399     CHECK("Failed to create region.");
00400     if (MB_SUCCESS != result) return result;
00401 
00402     result = mdbImpl->add_entities(this_gent, &new_handle, 1);
00403     CHECK("Failed to add region to geom set.");
00404     if (MB_SUCCESS != result) return result;
00405 
00406     if (file_id_tag)
00407       regions[i] = new_handle;
00408   
00409     n = fscanf(file_ptr, "%d ", &dum_int); CHECKN(1);
00410 
00411   } // end of reading regions
00412   
00413   if (file_id_tag) {
00414     result = readMeshIface->assign_ids( *file_id_tag, &regions[0], regions.size(), 1 );
00415     if (MB_SUCCESS != result)
00416       return result;
00417   }
00418 
00419   return MB_SUCCESS;
00420 }
00421 
00422 ErrorCode ReadSms::get_set(std::vector<EntityHandle> *sets,
00423                            int set_dim, int set_id,
00424                              Tag dim_tag,
00425                              EntityHandle &this_set,
00426                              const Tag* file_id_tag) 
00427 {
00428   ErrorCode result = MB_SUCCESS;
00429   
00430   if (set_dim < 0 || set_dim > 3)
00431     return MB_FILE_WRITE_ERROR;
00432   
00433   if ((int)sets[set_dim].size() <= set_id || 
00434       !sets[set_dim][set_id]) {
00435     if ((int)sets[set_dim].size() <= set_id) 
00436       sets[set_dim].resize(set_id+1, 0);
00437       
00438     if (!sets[set_dim][set_id]) {
00439       result = mdbImpl->create_meshset(MESHSET_SET, 
00440                                        sets[set_dim][set_id]);
00441       if (MB_SUCCESS != result) return result;
00442       result = mdbImpl->tag_set_data(globalId, 
00443                                      &sets[set_dim][set_id], 1,
00444                                      &set_id);
00445       if (MB_SUCCESS != result) return result;
00446       result = mdbImpl->tag_set_data(dim_tag,
00447                                      &sets[set_dim][set_id], 1,
00448                                      &set_dim);
00449       if (MB_SUCCESS != result) return result;
00450       
00451       if (file_id_tag) {
00452         result = mdbImpl->tag_set_data(*file_id_tag,
00453                                      &sets[set_dim][set_id], 1,
00454                                      &setId);
00455         ++setId;
00456       }
00457     }
00458   }
00459 
00460   this_set = sets[set_dim][set_id];
00461 
00462   return result;
00463 }
00464 
00465 ErrorCode ReadSms::read_parallel_info(FILE *file_ptr) 
00466 {
00467 //  ErrorCode result;
00468 
00469     // read partition info
00470   int nparts, part_id, num_ifaces, num_corner_ents;
00471   int num_read = fscanf(file_ptr, "%d %d %d %d", &nparts, &part_id, &num_ifaces, &num_corner_ents);
00472   if (!num_read) return MB_FAILURE;
00473   
00474     // read interfaces
00475   int iface_id, iface_dim, iface_own, num_iface_corners;
00476 //  EntityHandle this_iface;
00477   std::vector<int> *iface_corners;
00478   for (int i = 0; i < num_ifaces; i++) {
00479     num_read = fscanf(file_ptr, "%d %d %d %d", &iface_id, &iface_dim, &iface_own,
00480                       &num_iface_corners);
00481     if (!num_read) return MB_FAILURE;
00482     
00483 //    result = get_set(sets, iface_dim, iface_id, dim_tag, iface_own, this_iface);
00484 //    CHECK("Failed to make iface set.");
00485     
00486       // read the corner ids and store them on the set for now
00487     iface_corners = new std::vector<int>(num_iface_corners);
00488     for (int j = 0; j < num_iface_corners; j++) {
00489       num_read = fscanf(file_ptr, "%d", &(*iface_corners)[j]);
00490       if (!num_read) return MB_FAILURE;
00491     }
00492 
00493 //    result = tag_set_data(ifaceCornerTag, &this_iface, 1,
00494 //                          &iface_corners);
00495 //    CHECK("Failed to set iface corner tag.");
00496   }
00497 
00498     // interface data has been read
00499   return MB_SUCCESS;
00500 }
00501 
00502 ErrorCode ReadSms::add_entities( EntityHandle start,
00503                                    EntityHandle count,
00504                                    const Tag* file_id_tag )
00505 {
00506   if (!count || !file_id_tag)
00507     return MB_FAILURE;
00508 
00509   Range range;
00510   range.insert( start, start + count - 1 );
00511   return readMeshIface->assign_ids( *file_id_tag, range, 1 );
00512 }
00513 
00514 } // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines