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