 |
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 /**
2 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 * storing and accessing finite element mesh data.
4 *
5 * Copyright 2004 Sandia Corporation. Under the terms of Contract
6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 * retains certain rights in this software.
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 */
15
16 #ifdef WIN32
17 #ifdef _DEBUG
18 // turn off warnings that say they debugging identifier has been truncated
19 // this warning comes up when using some STL containers
20 #pragma warning( disable : 4786 )
21 #endif
22 #endif
23
24 #include "WriteAns.hpp"
25
26 #include <utility>
27 #include <algorithm>
28 #include <ctime>
29 #include <string>
30 #include <vector>
31 #include <cstdio>
32 #include <iostream>
33 #include <fstream>
34 #include <iomanip>
35
36 #include "moab/Interface.hpp"
37 #include "moab/Range.hpp"
38 #include <cassert>
39 #include "Internals.hpp"
40 #include "ExoIIUtil.hpp"
41 #include "MBTagConventions.hpp"
42
43 namespace moab
44 {
45
46 WriterIface* WriteAns::factory( Interface* iface )
47 {
48 return new WriteAns( iface );
49 }
50
51 WriteAns::WriteAns( Interface* impl ) : mbImpl( impl ), mCurrentMeshHandle( 0 ), mGlobalIdTag( 0 ), mMatSetIdTag( 0 )
52 {
53 assert( impl != NULL );
54
55 // impl->query_interface( mWriteIface );
56
57 // initialize in case tag_get_handle fails below
58 //! get and cache predefined tag handles
59 const int negone = -1;
60 impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
61 &negone );
62
63 impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
64 &negone );
65
66 impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
67 &negone );
68 }
69
70 WriteAns::~WriteAns()
71 {
72 // mbImpl->release_interface(mWriteIface);
73 }
74
75 ErrorCode WriteAns::write_file( const char* file_name,
76 const bool /* overwrite (commented out to remove warning) */,
77 const FileOptions&,
78 const EntityHandle* ent_handles,
79 const int num_sets,
80 const std::vector< std::string >&,
81 const Tag*,
82 int,
83 int )
84 {
85 assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
86
87 ErrorCode result;
88
89 // set SOLID45 element type to #60000, hope nobody has already...
90 const char* ETSolid45 = "60045";
91 const char* ETSolid92 = "60042";
92 const char* ETSolid95 = "60095";
93
94 // set Material id # to be used as default for all elements
95 // will need to be subsequently reassigned inside ANSYS
96 // Can, although have not, declare similar defaults for other attributes
97 const char* MATDefault = "1";
98
99 // create file streams for writing
100 std::ofstream node_file;
101 std::ofstream elem_file;
102 std::ofstream ans_file;
103
104 // get base filename from filename.ans
105 std::string temp_string;
106 std::string base_string;
107 base_string.assign( file_name );
108 base_string.replace( base_string.find_last_of( ".ans" ) - 3, 4, "" );
109
110 // open node file for writing
111 temp_string = base_string + ".node";
112 node_file.open( temp_string.c_str() );
113 node_file.setf( std::ios::scientific, std::ios::floatfield );
114 node_file.precision( 13 );
115
116 // open elem file for writing
117 temp_string = base_string + ".elem";
118 elem_file.open( temp_string.c_str() );
119
120 // open ans file for writing
121 ans_file.open( file_name );
122 ans_file << "/prep7" << std::endl;
123
124 // gather single output set
125 EntityHandle output_set = 0;
126 if( ent_handles && num_sets > 0 )
127 {
128 for( int i = 0; i < num_sets; i++ )
129 {
130 // from template, maybe can be removed
131 result = mbImpl->unite_meshset( output_set, ent_handles[i] );
132 if( result != MB_SUCCESS ) return result;
133 }
134 }
135
136 // search for all nodes
137 Range node_range;
138 result = mbImpl->get_entities_by_type( output_set, MBVERTEX, node_range, true );
139 if( result != MB_SUCCESS ) return result;
140
141 // Commented out until Seg Fault taken care of in gather_nodes...
142 // get any missing nodes which are needed for elements
143 // Range all_ent_range,missing_range;
144 // result=mbImpl->get_entities_by_handle(output_set,all_ent_range,true);
145 // if(result !=MB_SUCCESS) return result;
146 // result=mWriteIface->gather_nodes_from_elements(all_ent_range,0,missing_range);
147 // node_range.merge(missing_range);
148
149 // write the nodes
150 double coord[3];
151 for( Range::iterator it = node_range.begin(); it != node_range.end(); ++it )
152 {
153 EntityHandle node_handle = *it;
154
155 result = mbImpl->get_coords( &node_handle, 1, coord );
156 if( result != MB_SUCCESS ) return result;
157
158 node_file.width( 8 );
159 node_file << mbImpl->id_from_handle( node_handle );
160 node_file.width( 20 );
161 node_file << coord[0];
162 node_file.width( 20 );
163 node_file << coord[1];
164 node_file.width( 20 );
165 node_file << coord[2] << std::endl;
166 }
167
168 // update header to load nodes
169 ans_file << "nread," << base_string << ",node" << std::endl;
170
171 // search for all node sets (Dirichlet Sets)
172 Range node_mesh_sets;
173 int ns_id;
174 result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, node_mesh_sets );
175 if( result != MB_SUCCESS ) return result;
176
177 for( Range::iterator ns_it = node_mesh_sets.begin(); ns_it != node_mesh_sets.end(); ++ns_it )
178 {
179 result = mbImpl->tag_get_data( mDirichletSetTag, &( *ns_it ), 1, &ns_id );
180 if( result != MB_SUCCESS ) return result;
181 std::vector< EntityHandle > node_vector;
182 result = mbImpl->get_entities_by_handle( *ns_it, node_vector, true );
183 if( result != MB_SUCCESS ) return result;
184 // for every nodeset found, cycle through nodes in set:
185 for( std::vector< EntityHandle >::iterator node_it = node_vector.begin(); node_it != node_vector.end();
186 ++node_it )
187 {
188 int ns_node_id = mbImpl->id_from_handle( *node_it );
189 if( node_it == node_vector.begin() )
190 {
191 // select first node in new list
192 ans_file << "nsel,s,node,," << std::setw( 8 ) << ns_node_id << std::endl;
193 }
194 else
195 {
196 // append node to list
197 ans_file << "nsel,a,node,," << std::setw( 8 ) << ns_node_id << std::endl;
198 }
199 }
200 // create NS(#) nodeset
201 ans_file << "cm,NS" << ns_id << ",node" << std::endl;
202 }
203
204 // ANSYS Element format:
205 // I, J, K, L, M, N, O, P,etc... MAT, TYPE, REAL, SECNUM, ESYS, IEL
206 // I-P are nodes of element
207 // MAT = material number
208 // TYPE = Element type number
209 // REAL = Real constant set number
210 // SECNUM = section attribute number
211 // ESYS = coordinate system for nodes
212 // IEL = element # (unique?)
213 // For all nodes past 8, write on second line
214
215 // Write all MBTET elements
216 Range tet_range;
217 result = mbImpl->get_entities_by_type( output_set, MBTET, tet_range, true );
218 if( result != MB_SUCCESS ) return result;
219 for( Range::iterator elem_it = tet_range.begin(); elem_it != tet_range.end(); ++elem_it )
220 {
221 EntityHandle elem_handle = *elem_it;
222 int elem_id = mbImpl->id_from_handle( elem_handle );
223 std::vector< EntityHandle > conn;
224 result = mbImpl->get_connectivity( &elem_handle, 1, conn, false );
225 if( result != MB_SUCCESS ) return result;
226 // make sure 4 or 10 node tet
227 if( conn.size() != 4 && conn.size() != 10 )
228 {
229 std::cout << "Support not added for element type. \n";
230 return MB_FAILURE;
231 }
232 // write information for 4 node tet
233 if( conn.size() == 4 )
234 {
235 elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[1];
236 elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[2];
237 elem_file << std::setw( 8 ) << conn[3] << std::setw( 8 ) << conn[3];
238 elem_file << std::setw( 8 ) << conn[3] << std::setw( 8 ) << conn[3];
239
240 elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid45;
241 elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
242 elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
243 elem_file << std::endl;
244 }
245
246 // write information for 10 node tet
247 if( conn.size() == 10 )
248 {
249 elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[1];
250 elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[3];
251 elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[5];
252 elem_file << std::setw( 8 ) << conn[6] << std::setw( 8 ) << conn[7];
253
254 elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid92;
255 elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
256 elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
257 elem_file << std::endl;
258
259 elem_file << std::setw( 8 ) << conn[8] << std::setw( 8 ) << conn[9];
260 elem_file << std::endl;
261 }
262 }
263
264 // Write all MBHEX elements
265 Range hex_range;
266 result = mbImpl->get_entities_by_type( output_set, MBHEX, hex_range, true );
267 if( result != MB_SUCCESS ) return result;
268 for( Range::iterator elem_it = hex_range.begin(); elem_it != hex_range.end(); ++elem_it )
269 {
270 EntityHandle elem_handle = *elem_it;
271 int elem_id = mbImpl->id_from_handle( elem_handle );
272 std::vector< EntityHandle > conn;
273 result = mbImpl->get_connectivity( &elem_handle, 1, conn, false );
274 if( result != MB_SUCCESS ) return result;
275 // make sure supported hex type
276 if( conn.size() != 8 && conn.size() != 20 )
277 {
278 std::cout << "Support not added for element type. \n";
279 return MB_FAILURE;
280 }
281
282 // write information for 8 node hex
283 if( conn.size() == 8 )
284 {
285 elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[1];
286 elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[3];
287 elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[5];
288 elem_file << std::setw( 8 ) << conn[6] << std::setw( 8 ) << conn[7];
289
290 elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid45;
291 elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
292 elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
293 elem_file << std::endl;
294 }
295
296 // write information for 20 node hex
297 if( conn.size() == 20 )
298 {
299
300 elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[5];
301 elem_file << std::setw( 8 ) << conn[1] << std::setw( 8 ) << conn[0];
302 elem_file << std::setw( 8 ) << conn[7] << std::setw( 8 ) << conn[6];
303 elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[3];
304
305 elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid95;
306 elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
307 elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
308 elem_file << std::endl;
309
310 elem_file << std::setw( 8 ) << conn[16] << std::setw( 8 ) << conn[13];
311 elem_file << std::setw( 8 ) << conn[8] << std::setw( 8 ) << conn[12];
312 elem_file << std::setw( 8 ) << conn[18] << std::setw( 8 ) << conn[14];
313 elem_file << std::setw( 8 ) << conn[10] << std::setw( 8 ) << conn[15];
314 elem_file << std::setw( 8 ) << conn[19] << std::setw( 8 ) << conn[17];
315 elem_file << std::setw( 8 ) << conn[9] << std::setw( 8 ) << conn[11];
316 elem_file << std::endl;
317 }
318 }
319 // Write all MBPRISM elements
320 Range prism_range;
321 result = mbImpl->get_entities_by_type( output_set, MBPRISM, prism_range, true );
322 if( result != MB_SUCCESS ) return result;
323 for( Range::iterator elem_it = prism_range.begin(); elem_it != prism_range.end(); ++elem_it )
324 {
325 EntityHandle elem_handle = *elem_it;
326 int elem_id = mbImpl->id_from_handle( elem_handle );
327 std::vector< EntityHandle > conn;
328 result = mbImpl->get_connectivity( &elem_handle, 1, conn, false );
329 if( result != MB_SUCCESS ) return result;
330 // make sure supported prism type
331 if( conn.size() != 6 )
332 {
333 std::cout << "Support not added for element type. \n";
334 return MB_FAILURE;
335 }
336
337 // write information for 6 node prism
338 if( conn.size() == 6 )
339 {
340 elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[3];
341 elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[4];
342 elem_file << std::setw( 8 ) << conn[1] << std::setw( 8 ) << conn[2];
343 elem_file << std::setw( 8 ) << conn[5] << std::setw( 8 ) << conn[5];
344
345 elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid45;
346 elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
347 elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
348 elem_file << std::endl;
349 }
350 }
351
352 // create element types (for now writes all, even if not used)
353 ans_file << "et," << ETSolid45 << ",SOLID45" << std::endl;
354 ans_file << "et," << ETSolid92 << ",SOLID92" << std::endl;
355 ans_file << "et," << ETSolid95 << ",SOLID95" << std::endl;
356
357 // xxx pyramids, other elements later...
358
359 // write header to load elements
360 ans_file << "eread," << base_string << ",elem" << std::endl;
361
362 // search for all side sets (Neumann)
363 Range side_mesh_sets;
364 int ss_id;
365 result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, side_mesh_sets );
366 if( result != MB_SUCCESS ) return result;
367 // cycle through all sets found
368 for( Range::iterator ss_it = side_mesh_sets.begin(); ss_it != side_mesh_sets.end(); ++ss_it )
369 {
370 result = mbImpl->tag_get_data( mNeumannSetTag, &( *ss_it ), 1, &ss_id );
371 if( result != MB_SUCCESS ) return result;
372 std::vector< EntityHandle > elem_vector;
373 result = mbImpl->get_entities_by_handle( *ss_it, elem_vector, true );
374 if( result != MB_SUCCESS ) return result;
375
376 // cycle through elements in current side set
377 for( std::vector< EntityHandle >::iterator elem_it = elem_vector.begin(); elem_it != elem_vector.end();
378 ++elem_it )
379 {
380 EntityHandle elem_handle = *elem_it;
381
382 // instead of selecting current element in set, select its nodes...
383 std::vector< EntityHandle > conn;
384 result = mbImpl->get_connectivity( &elem_handle, 1, conn );
385 if( result != MB_SUCCESS ) return result;
386 if( elem_it == elem_vector.begin() )
387 {
388 ans_file << "nsel,s,node,," << std::setw( 8 ) << conn[0] << std::endl;
389 for( unsigned int i = 1; i < conn.size(); i++ )
390 {
391 ans_file << "nsel,a,node,," << std::setw( 8 ) << conn[i] << std::endl;
392 }
393 }
394 else
395 {
396 for( unsigned int i = 0; i < conn.size(); i++ )
397 {
398 ans_file << "nsel,a,node,," << std::setw( 8 ) << conn[i] << std::endl;
399 }
400 }
401 }
402 // create SS(#) node set
403 ans_file << "cm,SS" << ss_id << ",node" << std::endl;
404 }
405
406 // Gather all element blocks
407 Range matset;
408 int mat_id;
409 result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, matset );
410 if( result != MB_SUCCESS ) return result;
411 // cycle through all elem blocks
412 for( Range::iterator mat_it = matset.begin(); mat_it != matset.end(); ++mat_it )
413 {
414 EntityHandle matset_handle = *mat_it;
415 result = mbImpl->tag_get_data( mMaterialSetTag, &matset_handle, 1, &mat_id );
416 if( result != MB_SUCCESS ) return result;
417 std::vector< EntityHandle > mat_vector;
418 result = mbImpl->get_entities_by_handle( *mat_it, mat_vector, true );
419 if( result != MB_SUCCESS ) return result;
420 // cycle through elements in current mat set
421 for( std::vector< EntityHandle >::iterator elem_it = mat_vector.begin(); elem_it != mat_vector.end();
422 ++elem_it )
423 {
424 EntityHandle elem_handle = *elem_it;
425 int elem_id = mbImpl->id_from_handle( elem_handle );
426 if( elem_it == mat_vector.begin() )
427 {
428 ans_file << "esel,s,elem,," << std::setw( 8 ) << elem_id << std::endl;
429 }
430 else
431 {
432 ans_file << "esel,a,elem,," << std::setw( 8 ) << elem_id << std::endl;
433 }
434 }
435 // for each matset, write block command
436 ans_file << "cm,EB" << mat_id << ",elem" << std::endl;
437 }
438
439 // close all file streams
440 node_file.close();
441 elem_file.close();
442 ans_file.close();
443
444 return MB_SUCCESS;
445 }
446
447 } // namespace moab
1.9.1.