 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include <iostream>
2 #include <fstream>
3 #include <vector>
4 #include <cstdlib>
5 #include <sstream>
6 #include <cassert>
7
8 #include "ReadIDEAS.hpp"
9 #include "moab/Interface.hpp"
10 #include "Internals.hpp"
11 #include "moab/ReadUtilIface.hpp"
12 #include "FileTokenizer.hpp"
13 #include "MBTagConventions.hpp"
14 #include "moab/Range.hpp"
15 #include "moab/CN.hpp"
16
17 namespace moab
18 {
19
20 ReaderIface* ReadIDEAS::factory( Interface* iface )
21 {
22 return new ReadIDEAS( iface );
23 }
24
25 ReadIDEAS::ReadIDEAS( Interface* impl ) : MBI( impl )
26 {
27 impl->query_interface( readMeshIface );
28 }
29
30 ErrorCode ReadIDEAS::read_tag_values( const char* /* file_name */,
31 const char* /* tag_name */,
32 const FileOptions& /* opts */,
33 std::vector< int >& /* tag_values_out */,
34 const SubsetList* /* subset_list */ )
35 {
36 return MB_NOT_IMPLEMENTED;
37 }
38
39 ErrorCode ReadIDEAS::load_file( const char* fname,
40 const EntityHandle*,
41 const FileOptions& /*options*/,
42 const ReaderIface::SubsetList* subset_list,
43 const Tag* file_id_tag )
44 {
45 if( subset_list )
46 {
47 MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for IDEAS" );
48 }
49
50 file.open( fname );
51 if( !file.good() )
52 {
53 MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "Failed to open file: " << fname );
54 }
55
56 ErrorCode rval;
57
58 char line[10000];
59 file.getline( line, 10000 );
60 char* liter = line;
61 while( *liter && isspace( *liter ) )
62 ++liter;
63 if( *liter != '-' ) return MB_FAILURE;
64 ++liter;
65 if( *liter != '1' ) return MB_FAILURE;
66 while( *++liter )
67 if( !isspace( *liter ) ) return MB_FAILURE;
68
69 EntityHandle first_vertex = 0;
70 while( !file.eof() )
71 {
72 file.getline( line, 10000 );
73 unsigned int header_id = (unsigned int)strtol( line, NULL, 10 );
74
75 // Create vertices
76 if( DOUBLE_PRECISION_NODES0 == header_id || DOUBLE_PRECISION_NODES1 == header_id )
77 {
78 if( first_vertex ) // multiple vertex blocks?
79 return MB_FAILURE;
80 rval = create_vertices( first_vertex, file_id_tag );MB_CHK_SET_ERR( rval, "Failed to read vertices" );
81 }
82 // Create elements
83 else if( ELEMENTS0 == header_id || ELEMENTS1 == header_id || ELEMENTS2 == header_id )
84 {
85 if( !first_vertex ) // Need to read vertices first
86 return MB_FAILURE;
87 rval = create_elements( first_vertex, file_id_tag );MB_CHK_SET_ERR( rval, "Failed to read elements" );
88 }
89 // Skip everything else
90 else
91 {
92 rval = skip_header();
93 if( MB_SUCCESS != rval ) return MB_FAILURE;
94 }
95 }
96
97 file.close();
98 return MB_SUCCESS;
99 }
100
101 ErrorCode ReadIDEAS::skip_header()
102 {
103 // Go until finding a pair of -1 lines
104 char* ctmp;
105 char line[10000];
106 std::string s;
107
108 int end_of_block = 0;
109
110 long int il;
111
112 while( file.getline( line, 10000 ) )
113 {
114 il = std::strtol( line, &ctmp, 10 );
115 if( il == -1 )
116 {
117 s = ctmp;
118 if( s.empty() ) end_of_block++;
119 }
120 else
121 end_of_block = 0;
122
123 if( end_of_block >= 2 ) return MB_SUCCESS;
124 }
125
126 return MB_SUCCESS;
127 }
128
129 ErrorCode ReadIDEAS::create_vertices( EntityHandle& first_vertex, const Tag* file_id_tag )
130 {
131 // Read two lines: first has some data, second has coordinates
132 char line1[10000], line2[10000];
133 int il1, il2;
134 char *ctmp1, *ctmp2;
135 std::string s1, s2;
136
137 ErrorCode rval;
138
139 std::streampos top_of_block = file.tellg();
140 unsigned int num_verts = 0;
141
142 for( ;; )
143 {
144 // Read both lines
145 if( !file.getline( line1, 10000 ) ) return MB_FAILURE;
146 if( !file.getline( line2, 10000 ) ) return MB_FAILURE;
147
148 // Check if we are at the end of the block
149 il1 = std::strtol( line1, &ctmp1, 10 );
150 il2 = std::strtol( line2, &ctmp2, 10 );
151 if( ( il1 == -1 ) && ( il2 == -1 ) )
152 {
153 s1 = ctmp1;
154 s2 = ctmp2;
155 if( ( s1.empty() ) && ( s2.empty() ) ) break;
156 }
157 num_verts++;
158 }
159
160 file.seekg( top_of_block );
161
162 std::vector< double* > arrays;
163 rval = readMeshIface->get_node_coords( 3, num_verts, 0, first_vertex, arrays );
164 if( MB_SUCCESS != rval ) return rval;
165
166 Range verts;
167 verts.insert( first_vertex, first_vertex + num_verts - 1 );
168
169 double* x = arrays[0];
170 double* y = arrays[1];
171 double* z = arrays[2];
172
173 // For now, assume ids are sequential and begin with 1
174 Tag id_tag = MBI->globalId_tag();
175 const int beginning_node_id = 1;
176 int node_id = beginning_node_id;
177
178 for( unsigned int i = 0; i < num_verts; i++ )
179 {
180 if( !file.getline( line1, 10000 ) ) return MB_FAILURE;
181 if( !file.getline( line2, 10000 ) ) return MB_FAILURE;
182
183 // Get the id out of the 1st line. Check the assumption that node ids are
184 // sequential and begin with 1.
185 if( node_id != std::strtol( line1, &ctmp1, 10 ) )
186 MB_SET_ERR( MB_FAILURE, "node_id " << node_id << " line2:" << line2 << " ctmp1:" << ctmp1 );
187 else
188 ++node_id;
189
190 // Get the doubles out of the 2nd line
191 x[i] = std::strtod( line2, &ctmp2 );
192 y[i] = std::strtod( ctmp2 + 1, &ctmp2 );
193 z[i] = std::strtod( ctmp2 + 1, NULL );
194 }
195
196 if( !file.getline( line1, 10000 ) ) MB_SET_ERR( MB_FAILURE, " expect more lines" );
197 if( !file.getline( line2, 10000 ) ) MB_SET_ERR( MB_FAILURE, " expect more lines 2" );
198
199 // Tag the nodes with ids
200 rval = readMeshIface->assign_ids( id_tag, verts, beginning_node_id );MB_CHK_SET_ERR( rval, "Failed to assign IDs" );
201 if( file_id_tag )
202 {
203 rval = readMeshIface->assign_ids( *file_id_tag, verts, beginning_node_id );MB_CHK_SET_ERR( rval, "Failed to assign file IDs" );
204 }
205
206 return MB_SUCCESS;
207 }
208
209 ErrorCode ReadIDEAS::create_elements( EntityHandle vstart, const Tag* file_id_tag )
210 {
211 char line1[10000], line2[10000];
212 int il1, il2;
213 char *ctmp1, *ctmp2;
214 std::string s1, s2;
215 ErrorCode rval;
216 EntityHandle handle;
217
218 Tag mat_tag, phys_tag, id_tag;
219 rval = MBI->tag_get_handle( MAT_PROP_TABLE_TAG, 1, MB_TYPE_INTEGER, mat_tag, MB_TAG_DENSE | MB_TAG_CREAT );
220 if( MB_SUCCESS != rval && MB_ALREADY_ALLOCATED != rval ) return rval;
221 rval = MBI->tag_get_handle( PHYS_PROP_TABLE_TAG, 1, MB_TYPE_INTEGER, phys_tag, MB_TAG_DENSE | MB_TAG_CREAT );
222 if( MB_SUCCESS != rval && MB_ALREADY_ALLOCATED != rval ) return rval;
223 id_tag = MBI->globalId_tag();
224
225 for( ;; )
226 {
227 if( !file.getline( line1, 10000 ) || !file.getline( line2, 10000 ) ) return MB_FAILURE;
228
229 // Check if we are at the end of the block
230 il1 = std::strtol( line1, &ctmp1, 10 );
231 il2 = std::strtol( line2, &ctmp2, 10 );
232 if( ( il1 == -1 ) && ( il2 == -1 ) )
233 {
234 s1 = ctmp1;
235 s2 = ctmp2;
236 if( ( s1.empty() ) && ( s2.empty() ) ) return MB_SUCCESS;
237 }
238
239 // The first line describes attributes of the element other than connectivity.
240 const int element_id = strtol( line1 + 1, &ctmp1, 10 );
241 const int ideas_type = strtol( line1 + 11, &ctmp1, 10 );
242 const int phys_table = strtol( line1 + 21, &ctmp1, 10 );
243 const int mat_table = strtol( line1 + 31, &ctmp1, 10 );
244
245 // Determine the element type.
246 EntityType mb_type;
247 if( TRI0 == ideas_type || TRI1 == ideas_type )
248 mb_type = MBTRI;
249 else if( QUAD0 == ideas_type || QUAD1 == ideas_type )
250 mb_type = MBQUAD;
251 else if( TET == ideas_type )
252 mb_type = MBTET;
253 else if( HEX == ideas_type )
254 mb_type = MBHEX;
255 else if( WEDGE == ideas_type )
256 mb_type = MBPRISM;
257 else
258 {
259 std::cout << "IDEAS element type not yet added to MOAB reader." << std::endl;
260 return MB_NOT_IMPLEMENTED;
261 }
262
263 // Get the connectivity out of the 2nd line
264 std::stringstream ss( line2 );
265 const int n_conn = CN::VerticesPerEntity( mb_type );
266 EntityHandle conn[CN::MAX_NODES_PER_ELEMENT];
267 EntityHandle vert;
268 for( int i = 0; i < n_conn; ++i )
269 {
270 ss >> vert;
271 conn[i] = vstart + vert - 1;
272 }
273
274 // Make the element. According to the Gmsh 2.2.3 source code, the IDEAS
275 // canonical numbering is the same as MBCN.
276 rval = MBI->create_element( mb_type, conn, n_conn, handle );MB_CHK_SET_ERR( rval, "can't create elements of type " << mb_type );
277
278 // If the phys set does not already exist, create it.
279 Range phys_sets;
280 EntityHandle phys_set;
281 const void* const phys_set_id_val[] = { &phys_table };
282 rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &phys_tag, phys_set_id_val, 1, phys_sets );MB_CHK_SET_ERR( rval, "can't get phys sets" );
283 if( phys_sets.empty() )
284 {
285 rval = MBI->create_meshset( MESHSET_SET, phys_set );MB_CHK_SET_ERR( rval, "can't create phys set" );
286 rval = MBI->tag_set_data( phys_tag, &phys_set, 1, &phys_table );MB_CHK_SET_ERR( rval, "can't set tag to phys set" );
287 }
288 else if( 1 == phys_sets.size() )
289 {
290 phys_set = phys_sets.front();
291 }
292 else
293 {
294 return MB_MULTIPLE_ENTITIES_FOUND;
295 }
296 rval = MBI->add_entities( phys_set, &handle, 1 );MB_CHK_SET_ERR( rval, "can't add entities to phys set" );
297
298 // If the material set does not already exist, create it.
299 Range mat_sets;
300 EntityHandle mat_set;
301 const void* const mat_set_id_val[] = { &mat_table };
302 rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &mat_tag, mat_set_id_val, 1, mat_sets );
303 if( MB_SUCCESS != rval ) return rval;
304 if( mat_sets.empty() )
305 {
306 rval = MBI->create_meshset( MESHSET_SET, mat_set );
307 if( MB_SUCCESS != rval ) return rval;
308 rval = MBI->tag_set_data( mat_tag, &mat_set, 1, &mat_table );
309 if( MB_SUCCESS != rval ) return rval;
310 }
311 else if( 1 == mat_sets.size() )
312 {
313 mat_set = mat_sets.front();
314 }
315 else
316 {
317 return MB_MULTIPLE_ENTITIES_FOUND;
318 }
319 rval = MBI->add_entities( mat_set, &handle, 1 );
320 if( MB_SUCCESS != rval ) return rval;
321
322 // Tag the element with its id
323 rval = MBI->tag_set_data( id_tag, &handle, 1, &element_id );MB_CHK_SET_ERR( rval, "Failed to assign IDs" );
324 if( file_id_tag )
325 {
326 rval = MBI->tag_set_data( *file_id_tag, &handle, 1, &element_id );MB_CHK_SET_ERR( rval, "Failed to assign file IDs" );
327 }
328 }
329
330 return MB_SUCCESS;
331 }
332
333 } // namespace moab
1.9.1.