Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadUtil.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 "ReadUtil.hpp"
25 #include "moab/Core.hpp"
26 #include "AEntityFactory.hpp"
27 #include "moab/Error.hpp"
28 #include "SequenceManager.hpp"
29 #include "VertexSequence.hpp"
30 #include "ElementSequence.hpp"
31 
32 namespace moab
33 {
34 
35 #define RR \
36  if( MB_SUCCESS != result ) return result
37 
38 ReadUtil::ReadUtil( Core* mdb, Error* /*error_handler*/ ) : ReadUtilIface(), mMB( mdb ) {}
39 
40 ErrorCode ReadUtil::get_node_coords( const int /*num_arrays*/,
41  const int num_nodes,
42  const int preferred_start_id,
43  EntityHandle& actual_start_handle,
44  std::vector< double* >& arrays,
45  int sequence_size )
46 {
48  EntitySequence* seq = 0;
49 
50  if( num_nodes < 1 )
51  {
52  actual_start_handle = 0;
53  arrays.clear();
54  return MB_INDEX_OUT_OF_RANGE;
55  }
56 
57  // Create an entity sequence for these nodes
58  error = mMB->sequence_manager()->create_entity_sequence( MBVERTEX, num_nodes, 0, preferred_start_id,
59  actual_start_handle, seq, sequence_size );
60 
61  if( error != MB_SUCCESS ) return error;
62 
63  if( seq->start_handle() > actual_start_handle || seq->end_handle() < actual_start_handle ||
64  seq->end_handle() - actual_start_handle + 1 < (unsigned)num_nodes )
65  return MB_FAILURE;
66 
67  arrays.resize( 3 );
68 
69  error = static_cast< VertexSequence* >( seq )->get_coordinate_arrays( arrays[0], arrays[1], arrays[2] );
70  for( unsigned i = 0; i < arrays.size(); ++i )
71  if( arrays[i] ) arrays[i] += ( actual_start_handle - seq->start_handle() );
72 
73  return error;
74 }
75 
76 ErrorCode ReadUtil::get_element_connect( const int num_elements,
77  const int verts_per_element,
78  const EntityType mdb_type,
79  const int preferred_start_id,
80  EntityHandle& actual_start_handle,
81  EntityHandle*& array,
82  int sequence_size )
83 {
85  EntitySequence* seq;
86 
87  if( num_elements < 1 )
88  {
89  actual_start_handle = 0;
90  array = 0;
91  return MB_INDEX_OUT_OF_RANGE;
92  }
93 
94  // if (mdb_type <= MBVERTEX || mdb_type >= MBPOLYHEDRON || mdb_type == MBPOLYGON)
95  // return MB_TYPE_OUT_OF_RANGE;
96 
97  // Make an entity sequence to hold these elements
98  error =
99  mMB->sequence_manager()->create_entity_sequence( mdb_type, num_elements, verts_per_element, preferred_start_id,
100  actual_start_handle, seq, sequence_size );
101  if( MB_SUCCESS != error ) return error;
102 
103  if( seq->start_handle() > actual_start_handle || seq->end_handle() < actual_start_handle ||
104  seq->end_handle() - actual_start_handle + 1 < (unsigned)num_elements )
105  return MB_FAILURE;
106 
107  // Get an array for the connectivity
108  array = static_cast< ElementSequence* >( seq )->get_connectivity_array();
109  if( !array ) return MB_FAILURE;
110  array +=
111  ( actual_start_handle - seq->start_handle() ) * static_cast< ElementSequence* >( seq )->nodes_per_element();
112 
113  return error;
114 }
115 
117  const unsigned* flags,
118  EntityID start_id,
119  EntityHandle& start_handle )
120 {
121  if( num_sets < 1 )
122  {
123  start_handle = 0;
124  return MB_INDEX_OUT_OF_RANGE;
125  }
126 
128  EntitySequence* seq;
129  error = mMB->sequence_manager()->create_meshset_sequence( num_sets, start_id, flags, start_handle, seq );
130  if( MB_SUCCESS != error ) return error;
131 
132  if( seq->start_handle() > start_handle || seq->end_handle() < start_handle ||
133  seq->end_handle() - start_handle + 1 < (EntityHandle)num_sets )
134  return MB_FAILURE;
135 
136  return MB_SUCCESS;
137 }
138 
140  const int number_elements,
141  const int number_vertices_per_element,
142  const EntityHandle* conn_array )
143 {
144  EntityHandle tmp_hndl = start_handle;
145  AEntityFactory* adj_fact = mMB->a_entity_factory();
146 
147  // Iterate over the elements and update adjacency information
148  if( adj_fact != NULL && adj_fact->vert_elem_adjacencies() )
149  {
150  int j = 0;
151  for( int i = 0; i < number_elements; i++ )
152  {
153  adj_fact->notify_create_entity( tmp_hndl, ( conn_array + j ), number_vertices_per_element );
154  tmp_hndl++;
155  j += number_vertices_per_element;
156  }
157  }
158 
159  return MB_SUCCESS;
160 }
161 
162 ErrorCode ReadUtil::gather_related_ents( Range& partition, Range& related_ents, EntityHandle* file_set )
163 {
164  // Loop over any sets, getting contained ents
165  std::pair< Range::const_iterator, Range::const_iterator > pair_it = partition.equal_range( MBENTITYSET );
166 
167  ErrorCode result = MB_SUCCESS;
168  for( Range::const_iterator rit = pair_it.first; rit != pair_it.second; ++rit )
169  {
170  ErrorCode tmp_result = mMB->get_entities_by_handle( *rit, related_ents, Interface::UNION );
171  if( MB_SUCCESS != tmp_result ) result = tmp_result;
172  }
173  RR;
174 
175  // Gather adjacent ents of other dimensions
176  Range tmp_ents;
177  for( int dim = 3; dim >= 0; dim-- )
178  {
179  tmp_ents.clear();
180  ErrorCode tmp_result = mMB->get_adjacencies( related_ents, dim, false, tmp_ents, Interface::UNION );
181  if( MB_SUCCESS != tmp_result )
182  result = tmp_result;
183  else
184  related_ents.merge( tmp_ents );
185  }
186  RR;
187 
188  // Related ents includes the partition itself
189  related_ents.merge( partition );
190 
191  // Get contains-related sets
192  Range tmp_ents3, last_related;
193  if( file_set )
194  result = mMB->get_entities_by_type( *file_set, MBENTITYSET, tmp_ents3 );
195  else
196  result = mMB->get_entities_by_type( 0, MBENTITYSET, tmp_ents3 );RR;
197 
198  while( related_ents.size() != last_related.size() )
199  {
200  last_related = related_ents;
201  for( Range::iterator rit = tmp_ents3.begin(); rit != tmp_ents3.end(); ++rit )
202  {
203  if( related_ents.find( *rit ) != related_ents.end() ) continue;
204 
205  tmp_ents.clear();
206  result = mMB->get_entities_by_handle( *rit, tmp_ents, true );RR;
207  Range tmp_ents2 = intersect( tmp_ents, related_ents );
208 
209  // If the intersection is not empty, set is related
210  if( !tmp_ents2.empty() ) related_ents.insert( *rit );
211  }
212  }
213 
214  // Get child-related sets
215  last_related.clear();
216  while( related_ents.size() != last_related.size() )
217  {
218  last_related = related_ents;
219  std::pair< Range::const_iterator, Range::const_iterator > it_pair = last_related.equal_range( MBENTITYSET );
220 
221  for( Range::const_iterator rit = it_pair.first; rit != it_pair.second; ++rit )
222  {
223  // Get all children and add to related ents
224  tmp_ents.clear();
225  result = mMB->get_child_meshsets( *rit, tmp_ents, 0 );RR;
226  related_ents.merge( tmp_ents );
227  }
228  }
229 
230  // Get parent-related sets
231  last_related.clear();
232  while( related_ents.size() != last_related.size() )
233  {
234  last_related = related_ents;
235  std::pair< Range::const_iterator, Range::const_iterator > it_pair = last_related.equal_range( MBENTITYSET );
236 
237  for( Range::const_iterator rit = it_pair.first; rit != it_pair.second; ++rit )
238  {
239  // Get all parents and add to related ents
240  tmp_ents.clear();
241  result = mMB->get_parent_meshsets( *rit, tmp_ents, 0 );RR;
242  related_ents.merge( tmp_ents );
243  }
244  }
245 
246  return MB_SUCCESS;
247 }
248 
250  int* sense,
251  int bound_size,
252  int dim,
253  EntityHandle* bound_verts,
254  EntityType& etype )
255 {
256  // Get dimension of bounding entities
257  int bound_dim = CN::Dimension( TYPE_FROM_HANDLE( bound_ents[0] ) );
258  int indices[MAX_SUB_ENTITY_VERTICES];
259  const EntityHandle* connect = NULL;
260  std::vector< EntityHandle > tmp_connect;
261 
262  // Find the right entity type based on # bounding ents
263  int numv = 0, num_connect = 0;
264  ErrorCode result;
265  for( EntityType t = MBEDGE; t < MBENTITYSET; t++ )
266  {
267  int nindex = CN::NumSubEntities( t, bound_dim );
268  if( CN::Dimension( t ) != dim || nindex != bound_size ) continue;
269 
270  // Fill in vertices from bounding entity vertices
271  int nverts = CN::VerticesPerEntity( t );
272  std::fill( bound_verts, bound_verts + nverts, 0 );
273  for( int index = 0; index < nindex; index++ )
274  {
275  result = mMB->get_connectivity( bound_ents[index], connect, num_connect, false, &tmp_connect );
276  if( MB_SUCCESS != result ) return result;
277 
278  CN::SubEntityVertexIndices( t, bound_dim, index, indices );
279 
280  for( int c = 0; c < num_connect; c++ )
281  {
282  if( !bound_verts[indices[c]] )
283  {
284  bound_verts[indices[c]] = ( sense[index] > 0 ) ? connect[c] : connect[num_connect - c - 1];
285  numv++;
286  }
287  }
288  if( numv == nverts )
289  {
290  etype = t;
291  return MB_SUCCESS;
292  }
293  }
294  }
295 
296  // If we get here, we didn't get full connectivity
297  etype = MBMAXTYPE;
298  return MB_FAILURE;
299 }
300 
302 {
303  int size;
304  DataType type;
305  ErrorCode rval = mb->tag_get_bytes( tag, size );
306  if( MB_SUCCESS != rval ) return rval;
307  if( size != sizeof( int ) ) return MB_TYPE_OUT_OF_RANGE;
308  rval = mb->tag_get_data_type( tag, type );
309  if( type != MB_TYPE_OPAQUE && type != MB_TYPE_INTEGER ) return MB_TYPE_OUT_OF_RANGE;
310 
311  return MB_SUCCESS;
312 }
313 
314 ErrorCode ReadUtil::assign_ids( Tag id_tag, const Range& ents, int start )
315 {
316  ErrorCode rval = check_int_tag( mMB, id_tag );
317  if( MB_SUCCESS != rval ) return rval;
318 
319  Range tmp_range;
320  std::vector< int > data;
321  for( Range::const_pair_iterator i = ents.pair_begin(); i != ents.pair_end(); ++i )
322  {
323  data.resize( i->second + 1 - i->first );
324  for( std::vector< int >::iterator j = data.begin(); j != data.end(); ++j )
325  *j = start++;
326  tmp_range.clear();
327  tmp_range.insert( i->first, i->second );
328  rval = mMB->tag_set_data( id_tag, tmp_range, &data[0] );
329  if( MB_SUCCESS != rval ) return rval;
330  }
331 
332  return MB_SUCCESS;
333 }
334 
335 ErrorCode ReadUtil::assign_ids( Tag id_tag, const EntityHandle* ents, size_t num_ents, int start )
336 {
337  ErrorCode rval = check_int_tag( mMB, id_tag );
338  if( MB_SUCCESS != rval ) return rval;
339 
340  std::vector< int > data;
341  const EntityHandle* const end = ents + num_ents;
342  const EntityHandle* i = ents;
343  while( i != end )
344  {
345  const EntityHandle* next = std::find( i, end, 0u );
346  size_t size = next - i;
347  if( !size )
348  {
349  ++i;
350  continue;
351  }
352 
353  int id = start + ( i - ents );
354  data.resize( size );
355  for( std::vector< int >::iterator j = data.begin(); j != data.end(); ++j )
356  *j = id++;
357 
358  rval = mMB->tag_set_data( id_tag, i, size, &data[0] );
359  if( MB_SUCCESS != rval ) return rval;
360  }
361 
362  return MB_SUCCESS;
363 }
364 
366 {
367  ErrorCode rval = mMB->create_meshset( MESHSET_SET, gather_set );
368  if( MB_SUCCESS != rval ) return rval;
369 
370  Tag gather_set_tag;
371  rval = mMB->tag_get_handle( "GATHER_SET", 1, MB_TYPE_INTEGER, gather_set_tag, MB_TAG_CREAT | MB_TAG_SPARSE );
372  if( MB_SUCCESS != rval ) return rval;
373 
374  int gather_val = 1;
375  rval = mMB->tag_set_data( gather_set_tag, &gather_set, 1, &gather_val );
376  if( MB_SUCCESS != rval ) return rval;
377 
378  return MB_SUCCESS;
379 }
380 
382 {
383  Tag gather_set_tag;
384  ErrorCode rval = mMB->tag_get_handle( "GATHER_SET", 1, MB_TYPE_INTEGER, gather_set_tag, MB_TAG_SPARSE );
385  if( MB_SUCCESS != rval ) return rval;
386 
387  int gather_val = 1;
388  void* vals[] = { &gather_val };
389  Range gather_sets;
390  rval = mMB->get_entities_by_type_and_tag( 0, MBENTITYSET, &gather_set_tag, vals, 1, gather_sets );
391  if( MB_SUCCESS != rval ) return rval;
392 
393  if( gather_sets.empty() ) return MB_ENTITY_NOT_FOUND;
394 
395  gather_set = gather_sets[0];
396 
397  return MB_SUCCESS;
398 }
399 
400 } // namespace moab