Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
WriteAns.cpp
Go to the documentation of this file.
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