moab
|
00001 #include "WriteGmsh.hpp" 00002 #include "moab/CN.hpp" 00003 #include "MBTagConventions.hpp" 00004 #include "MBParallelConventions.h" 00005 #include "moab/Interface.hpp" 00006 #include "moab/Range.hpp" 00007 #include "moab/WriteUtilIface.hpp" 00008 #include "moab/FileOptions.hpp" 00009 #include "GmshUtil.hpp" 00010 00011 #include <fstream> 00012 #include <map> 00013 #include <set> 00014 00015 namespace moab { 00016 00017 const int DEFAULT_PRECISION = 10; 00018 00019 WriterIface *WriteGmsh::factory( Interface* iface ) 00020 { return new WriteGmsh( iface ); } 00021 00022 WriteGmsh::WriteGmsh(Interface *impl) 00023 : mbImpl(impl) 00024 { 00025 impl->query_interface(mWriteIface); 00026 } 00027 00028 WriteGmsh::~WriteGmsh() 00029 { 00030 mbImpl->release_interface(mWriteIface); 00031 } 00032 00033 00034 // A structure to store per-element information. 00035 struct ElemInfo { 00036 void set( int st, int idt ) { 00037 while (count < st) 00038 sets[count++] = 0; 00039 sets[count++] = idt; 00040 } 00041 int count; // number of valid entries in sets[] 00042 int sets[3]; // ids of owning block, geom, and partition; respectively 00043 int id; // global ID of element 00044 int type; // Gmsh element type 00045 }; 00046 00048 ErrorCode WriteGmsh::write_file(const char *file_name, 00049 const bool overwrite, 00050 const FileOptions& options, 00051 const EntityHandle *output_list, 00052 const int num_sets, 00053 const std::vector<std::string>& , 00054 const Tag* , 00055 int , 00056 int ) 00057 { 00058 ErrorCode rval; 00059 Tag global_id = 0, block_tag = 0, geom_tag = 0, prtn_tag = 0; 00060 00061 if (!overwrite) 00062 { 00063 rval = mWriteIface->check_doesnt_exist( file_name ); 00064 if (MB_SUCCESS != rval) 00065 return rval; 00066 } 00067 00068 // Get tags 00069 mbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, global_id ); 00070 mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, block_tag ); 00071 if (global_id) 00072 mbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag ); 00073 mbImpl->tag_get_handle( PARALLEL_PARTITION_TAG_NAME, 1, MB_TYPE_INTEGER, prtn_tag ); 00074 00075 00076 // Define arrays to hold entity sets of interest 00077 Range sets[3]; 00078 Tag set_tags[] = { block_tag, geom_tag, prtn_tag }; 00079 Tag set_ids[] = { block_tag, 0 /*global_id*/, prtn_tag }; 00080 00081 // Get entities to write 00082 Range elements, nodes; 00083 if (!output_list) 00084 { 00085 rval = mbImpl->get_entities_by_dimension( 0, 0, nodes, false ); 00086 if (MB_SUCCESS != rval) 00087 return rval; 00088 for (int d = 1; d < 3; ++d) 00089 { 00090 Range tmp_range; 00091 rval = mbImpl->get_entities_by_dimension( 0, d, tmp_range, false ); 00092 if (MB_SUCCESS != rval) 00093 return rval; 00094 elements.merge( tmp_range ); 00095 } 00096 00097 for (int s = 0; s < 3; ++s) 00098 if (set_tags[s]) 00099 { 00100 rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, set_tags+s, 0, 1, sets[s] ); 00101 if (MB_SUCCESS != rval) return rval; 00102 } 00103 } 00104 else 00105 { 00106 for (int i = 0; i < num_sets; ++i) 00107 { 00108 EntityHandle set = output_list[i]; 00109 for (int d = 1; d < 3; ++d) 00110 { 00111 Range tmp_range, tmp_nodes; 00112 rval = mbImpl->get_entities_by_dimension( set, d, tmp_range, true ); 00113 if (rval != MB_SUCCESS) 00114 return rval; 00115 elements.merge( tmp_range ); 00116 rval = mbImpl->get_adjacencies( tmp_range, set, false, tmp_nodes ); 00117 if (rval != MB_SUCCESS) 00118 return rval; 00119 nodes.merge( tmp_nodes ); 00120 } 00121 00122 for (int s = 0; s < 3; ++s) 00123 if (set_tags[s]) 00124 { 00125 Range tmp_range; 00126 rval = mbImpl->get_entities_by_type_and_tag( set, MBENTITYSET, set_tags+s, 0, 1, tmp_range ); 00127 if (MB_SUCCESS != rval) return rval; 00128 sets[s].merge( tmp_range ); 00129 int junk; 00130 rval = mbImpl->tag_get_data( set_tags[s], &set, 1, &junk ); 00131 if (MB_SUCCESS == rval) 00132 sets[s].insert( set ); 00133 } 00134 } 00135 } 00136 00137 if (elements.empty()) 00138 { 00139 mWriteIface->report_error( "Nothing to write.\n" ); 00140 return MB_ENTITY_NOT_FOUND; 00141 } 00142 00143 // get global IDs for all elements. 00144 // First try to get from tag. If tag is not defined or not set 00145 // for all elements, use handle value instead. 00146 std::vector<int> global_id_array(elements.size()); 00147 std::vector<int>::iterator id_iter; 00148 if (!global_id || MB_SUCCESS != 00149 mbImpl->tag_get_data( global_id, elements, &global_id_array[0] ) ) 00150 { 00151 id_iter = global_id_array.begin(); 00152 for (Range::iterator i = elements.begin(); i != elements.end(); ++i, ++id_iter) 00153 *id_iter = mbImpl->id_from_handle( *i ); 00154 } 00155 00156 // Figure out the maximum ID value so we know where to start allocating 00157 // new IDs when we encounter ID conflits. 00158 int max_id = 0; 00159 for (id_iter = global_id_array.begin(); id_iter != global_id_array.end(); ++id_iter) 00160 if (*id_iter > max_id) 00161 max_id = *id_iter; 00162 00163 // Initialize ElemInfo struct for each element 00164 std::map<EntityHandle,ElemInfo> elem_sets; // per-element info 00165 std::set<int> elem_global_ids; // temporary for finding duplicate IDs 00166 id_iter = global_id_array.begin(); 00167 // Iterate backwards to give highest-dimension entities first dibs for 00168 // a conflicting ID. 00169 for (Range::reverse_iterator i = elements.rbegin(); i != elements.rend(); ++i) 00170 { 00171 int id = *id_iter; ++id_iter; 00172 if (!elem_global_ids.insert(id).second) 00173 id = max_id++; 00174 00175 ElemInfo& ei = elem_sets[*i]; 00176 ei.count = 0; 00177 ei.id = id; 00178 00179 EntityType type = mbImpl->type_from_handle( *i ); 00180 int num_vtx; 00181 const EntityHandle* conn; 00182 rval = mbImpl->get_connectivity( *i, conn, num_vtx ); 00183 if (MB_SUCCESS != rval) 00184 return rval; 00185 00186 ei.type = GmshUtil::get_gmsh_type( type, num_vtx ); 00187 if (ei.type < 0) 00188 { 00189 mWriteIface->report_error( "Gmem file format does not support element " 00190 " of type %s with %d vertices.\n", 00191 CN::EntityTypeName( type ), num_vtx ); 00192 return MB_FILE_WRITE_ERROR; 00193 } 00194 } 00195 // Don't need these any more, free memory. 00196 elem_global_ids.clear(); 00197 global_id_array.clear(); 00198 00199 // For each material set, geometry set, or partition; store 00200 // the ID of the set on each element. 00201 for (int s = 0; s < 3; ++s) 00202 { 00203 if (!set_tags[s]) 00204 continue; 00205 00206 for (Range::iterator i = sets[s].begin(); i != sets[s].end(); ++i) 00207 { 00208 int id; 00209 if (set_ids[s]) 00210 { 00211 rval = mbImpl->tag_get_data( set_ids[s], &*i, 1, &id ); 00212 if (MB_SUCCESS != rval) 00213 return rval; 00214 } 00215 else 00216 { 00217 id = mbImpl->id_from_handle( *i ); 00218 } 00219 00220 Range elems; 00221 rval = mbImpl->get_entities_by_handle( *i, elems ); 00222 if (MB_SUCCESS != rval) 00223 return rval; 00224 00225 elems = intersect( elems, elements ); 00226 for (Range::iterator j = elems.begin(); j != elems.end(); ++j) 00227 elem_sets[*j].set( s, id ); 00228 } 00229 } 00230 00231 00232 // Create file 00233 std::ofstream out( file_name ); 00234 if (!out) 00235 return MB_FILE_DOES_NOT_EXIST; 00236 00237 00238 // Write header 00239 out << "$MeshFormat" << std::endl; 00240 out << "2.0 0 " << sizeof(double) << std::endl; 00241 out << "$EndMeshFormat" << std::endl; 00242 00243 // Set precision for node coordinates 00244 int precision; 00245 if (MB_SUCCESS != options.get_int_option( "PRECISION", precision )) 00246 precision = DEFAULT_PRECISION; 00247 const int old_precision = out.precision(); 00248 out.precision( precision ); 00249 00250 // Write nodes 00251 out << "$Nodes" << std::endl; 00252 out << nodes.size() << std::endl; 00253 std::vector<double> coords(3*nodes.size()); 00254 rval = mbImpl->get_coords( nodes, &coords[0] ); 00255 if (MB_SUCCESS != rval) 00256 return rval; 00257 std::vector<double>::iterator c = coords.begin(); 00258 for (Range::iterator i = nodes.begin(); i != nodes.end(); ++i) 00259 { 00260 out << mbImpl->id_from_handle( *i ); 00261 out << " " << *c; ++c; 00262 out << " " << *c; ++c; 00263 out << " " << *c; ++c; 00264 out << std::endl; 00265 } 00266 out << "$EndNodes" << std::endl; 00267 coords.clear(); 00268 00269 // Restore stream state 00270 out.precision( old_precision ); 00271 00272 // Write elements 00273 out << "$Elements" << std::endl; 00274 out << elem_sets.size() << std::endl; 00275 for (std::map<EntityHandle,ElemInfo>::iterator i = elem_sets.begin(); 00276 i != elem_sets.end(); ++i) 00277 { 00278 int num_vtx; 00279 const EntityHandle* conn; 00280 rval = mbImpl->get_connectivity( i->first, conn, num_vtx ); 00281 if (MB_SUCCESS != rval) 00282 return rval; 00283 out << i->second.id << ' ' << i->second.type << ' ' << i->second.count; 00284 for (int j = 0; j < i->second.count; ++j) 00285 out << ' ' << i->second.sets[j]; 00286 00287 const int* order = GmshUtil::gmshElemTypes[i->second.type].node_order; 00288 00289 // need to re-order vertices 00290 if (order) 00291 { 00292 for (int j = 0; j < num_vtx; ++j) 00293 out << ' ' << mbImpl->id_from_handle( conn[order[j]] ); 00294 } 00295 else 00296 { 00297 for (int j = 0; j < num_vtx; ++j) 00298 out << ' ' << mbImpl->id_from_handle( conn[j] ); 00299 } 00300 out << std::endl; 00301 } 00302 out << "$EndElements" << std::endl; 00303 00304 00305 // done 00306 return MB_SUCCESS; 00307 } 00308 00309 } // namespace moab 00310 00311