Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
HigherOrderFactory.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 "moab/HigherOrderFactory.hpp"
25 #include "SequenceManager.hpp"
26 #include "UnstructuredElemSeq.hpp"
27 #include "VertexSequence.hpp"
28 #include "AEntityFactory.hpp"
29 #include "moab/Core.hpp"
30 #include "moab/CN.hpp"
31 #include <cassert>
32 #include <algorithm>
33 
34 namespace moab
35 {
36 
37 using namespace std;
38 
39 HigherOrderFactory::HigherOrderFactory( Core* MB, Interface::HONodeAddedRemoved* function_object )
40  : mMB( MB ), mHONodeAddedRemoved( function_object )
41 {
42  initialize_map();
43 }
44 HigherOrderFactory::~HigherOrderFactory() {}
45 
46 // bool HigherOrderFactory::mMapInitialized = false;
47 
48 void HigherOrderFactory::initialize_map()
49 {
50  // if(mMapInitialized)
51  // return;
52 
53  for( EntityType i = MBVERTEX; i < MBMAXTYPE; i++ )
54  {
55  const CN::ConnMap& canon_map = CN::mConnectivityMap[i][0];
56  unsigned char( &this_map )[8][8] = mNodeMap[i];
57  int num_node = CN::VerticesPerEntity( i );
58  for( int j = 0; j < canon_map.num_sub_elements; j++ )
59  {
60  unsigned char x = canon_map.conn[j][0];
61  unsigned char y = canon_map.conn[j][1];
62  this_map[x][y] = num_node;
63  this_map[y][x] = num_node;
64  num_node++;
65  }
66  }
67 
68  // mMapInitialized = true;
69 }
70 
71 ErrorCode HigherOrderFactory::convert( const EntityHandle meshset,
72  const bool mid_edge_nodes,
73  const bool mid_face_nodes,
74  const bool mid_volume_nodes )
75 {
76  Range entities;
77  mMB->get_entities_by_handle( meshset, entities, true );
78  return convert( entities, mid_edge_nodes, mid_face_nodes, mid_volume_nodes );
79 }
80 
81 ErrorCode HigherOrderFactory::convert( const Range& entities,
82  const bool mid_edge_nodes,
83  const bool mid_face_nodes,
84  const bool mid_volume_nodes )
85 
86 {
87 
88  // TODO -- add some more code to prevent from splitting of entity sequences when we don't need
89  // to. Say we have all hex8's in our mesh and 3 falses are passed in. In the end, no conversion
90  // will happen, but the sequences could still be split up.
91 
92  // find out what entity sequences we need to convert
93  // and where, if necessary, to split them
94 
95  SequenceManager* seq_manager = mMB->sequence_manager();
96  Range::const_pair_iterator p_iter;
97  for( p_iter = entities.const_pair_begin(); p_iter != entities.const_pair_end(); ++p_iter )
98  {
99 
100  EntityHandle h = p_iter->first;
101  while( h <= p_iter->second )
102  {
103 
104  EntitySequence* seq;
105  ErrorCode rval = seq_manager->find( h, seq );
106  if( MB_SUCCESS != rval ) return rval;
107 
108  if( seq->type() == MBVERTEX || seq->type() >= MBENTITYSET ) return MB_TYPE_OUT_OF_RANGE;
109 
110  // make sequence is not structured mesh
111  ElementSequence* elemseq = static_cast< ElementSequence* >( seq );
112  if( NULL == elemseq->get_connectivity_array() ) return MB_NOT_IMPLEMENTED;
113 
114  EntityHandle last = p_iter->second;
115  if( last > seq->end_handle() ) last = seq->end_handle();
116 
117  rval = convert_sequence( elemseq, h, last, mid_edge_nodes, mid_face_nodes, mid_volume_nodes );
118  if( MB_SUCCESS != rval ) return rval;
119 
120  h = last + 1;
121  }
122  }
123 
124  return MB_SUCCESS;
125 }
126 
127 ErrorCode HigherOrderFactory::convert_sequence( ElementSequence* seq,
128  EntityHandle start,
129  EntityHandle end,
130  bool mid_edge_nodes,
131  bool mid_face_nodes,
132  bool mid_volume_nodes )
133 {
134 
135  ErrorCode status = MB_SUCCESS;
136 
137  // lets make sure parameters are ok before we continue
138  switch( seq->type() )
139  {
140  default:
141  return MB_TYPE_OUT_OF_RANGE;
142  case MBEDGE:
143  mid_face_nodes = false;
144  mid_volume_nodes = false;
145  break;
146  case MBTRI:
147  case MBQUAD:
148  mid_volume_nodes = false;
149  break;
150  case MBTET:
151  case MBHEX:
152  case MBPRISM:
153  case MBPYRAMID:
154  case MBKNIFE:
155  break;
156  }
157 
158  // calculate number of nodes in target configuration
159  unsigned nodes_per_elem = CN::VerticesPerEntity( seq->type() );
160  if( mid_edge_nodes ) nodes_per_elem += ( seq->type() == MBEDGE ) ? 1 : CN::NumSubEntities( seq->type(), 1 );
161  if( mid_face_nodes )
162  nodes_per_elem += ( CN::Dimension( seq->type() ) == 2 ) ? 1 : CN::NumSubEntities( seq->type(), 2 );
163  if( mid_volume_nodes ) nodes_per_elem += 1;
164 
165  if( nodes_per_elem == seq->nodes_per_element() ) return MB_SUCCESS;
166 
167  Tag deletable_nodes;
168  status = mMB->tag_get_handle( 0, 1, MB_TYPE_BIT, deletable_nodes, MB_TAG_CREAT | MB_TAG_BIT );
169  if( MB_SUCCESS != status ) return status;
170 
171  UnstructuredElemSeq* new_seq = new UnstructuredElemSeq( start, end - start + 1, nodes_per_elem, end - start + 1 );
172 
173  copy_corner_nodes( seq, new_seq );
174 
175  if( seq->has_mid_edge_nodes() && mid_edge_nodes )
176  status = copy_mid_edge_nodes( seq, new_seq );
177  else if( seq->has_mid_edge_nodes() && !mid_edge_nodes )
178  status = remove_mid_edge_nodes( seq, start, end, deletable_nodes );
179  else if( !seq->has_mid_edge_nodes() && mid_edge_nodes )
180  status = zero_mid_edge_nodes( new_seq );
181  if( MB_SUCCESS != status ) return status;
182 
183  if( seq->has_mid_face_nodes() && mid_face_nodes )
184  status = copy_mid_face_nodes( seq, new_seq );
185  else if( seq->has_mid_face_nodes() && !mid_face_nodes )
186  status = remove_mid_face_nodes( seq, start, end, deletable_nodes );
187  else if( !seq->has_mid_face_nodes() && mid_face_nodes )
188  status = zero_mid_face_nodes( new_seq );
189  if( MB_SUCCESS != status )
190  {
191  mMB->tag_delete( deletable_nodes );
192  return status;
193  }
194 
195  if( seq->has_mid_volume_nodes() && mid_volume_nodes )
196  status = copy_mid_volume_nodes( seq, new_seq );
197  else if( seq->has_mid_volume_nodes() && !mid_volume_nodes )
198  status = remove_mid_volume_nodes( seq, start, end, deletable_nodes );
199  else if( !seq->has_mid_volume_nodes() && mid_volume_nodes )
200  status = zero_mid_volume_nodes( new_seq );
201  if( MB_SUCCESS != status )
202  {
203  mMB->tag_delete( deletable_nodes );
204  return status;
205  }
206 
207  // gather nodes that were marked
208  Range nodes;
209  mMB->get_entities_by_type_and_tag( 0, MBVERTEX, &deletable_nodes, NULL, 1, nodes );
210 
211  // EntityHandle low_meshset;
212  // int dum;
213  // low_meshset = CREATE_HANDLE(MBENTITYSET, 0, dum);
214 
215  for( Range::iterator iter = nodes.begin(); iter != nodes.end(); ++iter )
216  {
217  unsigned char marked = 0;
218  mMB->tag_get_data( deletable_nodes, &( *iter ), 1, &marked );
219  if( marked )
220  {
221  // we can delete it
222  if( mHONodeAddedRemoved ) mHONodeAddedRemoved->node_removed( *iter );
223  mMB->delete_entities( &( *iter ), 1 );
224  }
225  }
226 
227  const bool create_midedge = !seq->has_mid_edge_nodes() && mid_edge_nodes;
228  const bool create_midface = !seq->has_mid_face_nodes() && mid_face_nodes;
229  const bool create_midvolm = !seq->has_mid_volume_nodes() && mid_volume_nodes;
230 
231  mMB->tag_delete( deletable_nodes );
232 
233  status = mMB->sequence_manager()->replace_subsequence( new_seq );
234  if( MB_SUCCESS != status )
235  {
236  SequenceData* data = new_seq->data();
237  delete new_seq;
238  delete data;
239  return status;
240  }
241 
242  if( create_midedge )
243  {
244  status = add_mid_edge_nodes( new_seq );
245  if( MB_SUCCESS != status ) return status;
246  }
247  if( create_midface )
248  {
249  status = add_mid_face_nodes( new_seq );
250  if( MB_SUCCESS != status ) return status;
251  }
252  if( create_midvolm )
253  {
254  status = add_mid_volume_nodes( new_seq );
255  if( MB_SUCCESS != status ) return status;
256  }
257 
258  return status;
259 }
260 
261 ErrorCode HigherOrderFactory::add_mid_volume_nodes( ElementSequence* seq )
262 {
263  EntityType this_type = seq->type();
264  SequenceManager* seq_manager = mMB->sequence_manager();
265 
266  // find out where in the connectivity list to add these new mid volume nodes
267  int edge_factor = seq->has_mid_edge_nodes() ? 1 : 0;
268  int face_factor = seq->has_mid_face_nodes() ? 1 : 0;
269  // offset by number of higher order nodes on edges if they exist
270  int num_corner_nodes = CN::VerticesPerEntity( this_type );
271  int new_node_index = num_corner_nodes;
272  new_node_index += edge_factor * CN::mConnectivityMap[this_type][0].num_sub_elements;
273  new_node_index += face_factor * CN::mConnectivityMap[this_type][1].num_sub_elements;
274 
275  EntityHandle* element = seq->get_connectivity_array();
276  EntityHandle curr_handle = seq->start_handle();
277  int nodes_per_element = seq->nodes_per_element();
278  EntityHandle* end_element = element + nodes_per_element * ( seq->size() );
279 
280  // iterate over the elements
281  for( ; element < end_element; element += nodes_per_element )
282  {
283  // find the centroid of this element
284  double tmp_coords[3], sum_coords[3] = { 0, 0, 0 };
285  EntitySequence* eseq = NULL;
286  for( int i = 0; i < num_corner_nodes; i++ )
287  {
288  seq_manager->find( element[i], eseq );
289  static_cast< VertexSequence* >( eseq )->get_coordinates( element[i], tmp_coords[0], tmp_coords[1],
290  tmp_coords[2] );
291  sum_coords[0] += tmp_coords[0];
292  sum_coords[1] += tmp_coords[1];
293  sum_coords[2] += tmp_coords[2];
294  }
295  sum_coords[0] /= num_corner_nodes;
296  sum_coords[1] /= num_corner_nodes;
297  sum_coords[2] /= num_corner_nodes;
298 
299  // create a new vertex at the centroid
300  mMB->create_vertex( sum_coords, element[new_node_index] );
301 
302  if( mHONodeAddedRemoved ) mHONodeAddedRemoved->node_added( element[new_node_index], curr_handle );
303 
304  curr_handle++;
305  }
306 
307  return MB_SUCCESS;
308 }
309 
310 ErrorCode HigherOrderFactory::add_mid_face_nodes( ElementSequence* seq )
311 {
312  EntityType this_type = seq->type();
313  SequenceManager* seq_manager = mMB->sequence_manager();
314  int num_vertices = CN::VerticesPerEntity( this_type );
315  int num_edges = CN::mConnectivityMap[this_type][0].num_sub_elements;
316  num_edges = seq->has_mid_edge_nodes() ? num_edges : 0;
317  int num_faces = CN::mConnectivityMap[this_type][1].num_sub_elements;
318 
319  const CN::ConnMap& entity_faces = CN::mConnectivityMap[this_type][1];
320 
321  EntityHandle* element = seq->get_connectivity_array();
322  EntityHandle curr_handle = seq->start_handle();
323  int nodes_per_element = seq->nodes_per_element();
324  EntityHandle* end_element = element + nodes_per_element * ( seq->size() );
325 
326  EntityHandle tmp_face_conn[4]; // max face nodes = 4
327  std::vector< EntityHandle > adjacent_entities( 4 );
328 
329  double tmp_coords[3];
330 
331  // iterate over the elements
332  for( ; element < end_element; element += nodes_per_element )
333  {
334  // for each edge in this entity
335  for( int i = 0; i < num_faces; i++ )
336  {
337  // a node was already assigned
338  if( element[i + num_edges + num_vertices] != 0 ) continue;
339 
340  tmp_face_conn[0] = element[entity_faces.conn[i][0]];
341  tmp_face_conn[1] = element[entity_faces.conn[i][1]];
342  tmp_face_conn[2] = element[entity_faces.conn[i][2]];
343  if( entity_faces.num_corners_per_sub_element[i] == 4 )
344  tmp_face_conn[3] = element[entity_faces.conn[i][3]];
345  else
346  tmp_face_conn[3] = 0;
347 
348  EntityHandle already_made_node = center_node_exist( tmp_face_conn, adjacent_entities );
349 
350  if( already_made_node )
351  {
352  element[i + num_edges + num_vertices] = already_made_node;
353  }
354  // create a node
355  else
356  {
357  EntitySequence* tmp_sequence = NULL;
358  double sum_coords[3] = { 0, 0, 0 };
359  int max_nodes = entity_faces.num_corners_per_sub_element[i];
360  for( int k = 0; k < max_nodes; k++ )
361  {
362  seq_manager->find( tmp_face_conn[k], tmp_sequence );
363  static_cast< VertexSequence* >( tmp_sequence )
364  ->get_coordinates( tmp_face_conn[k], tmp_coords[0], tmp_coords[1], tmp_coords[2] );
365  sum_coords[0] += tmp_coords[0];
366  sum_coords[1] += tmp_coords[1];
367  sum_coords[2] += tmp_coords[2];
368  }
369 
370  sum_coords[0] /= max_nodes;
371  sum_coords[1] /= max_nodes;
372  sum_coords[2] /= max_nodes;
373 
374  mMB->create_vertex( sum_coords, element[i + num_edges + num_vertices] );
375  }
376 
377  if( mHONodeAddedRemoved )
378  mHONodeAddedRemoved->node_added( element[i + num_edges + num_vertices], curr_handle );
379  }
380 
381  curr_handle++;
382  }
383 
384  return MB_SUCCESS;
385 }
386 
387 ErrorCode HigherOrderFactory::add_mid_edge_nodes( ElementSequence* seq )
388 {
389  // for each node, need to see if it was already created.
390  EntityType this_type = seq->type();
391  SequenceManager* seq_manager = mMB->sequence_manager();
392 
393  // offset by number of corner nodes
394  int num_vertices = CN::VerticesPerEntity( this_type );
395  int num_edges = CN::mConnectivityMap[this_type][0].num_sub_elements;
396 
397  const CN::ConnMap& entity_edges = CN::mConnectivityMap[this_type][0];
398 
399  EntityHandle* element = seq->get_connectivity_array();
400  EntityHandle curr_handle = seq->start_handle();
401  int nodes_per_element = seq->nodes_per_element();
402  EntityHandle* end_element = element + nodes_per_element * ( seq->size() );
403 
404  EntityHandle tmp_edge_conn[2];
405  std::vector< EntityHandle > adjacent_entities( 32 );
406 
407  double tmp_coords[3];
408 
409  // iterate over the elements
410  for( ; element < end_element; element += nodes_per_element )
411  {
412  // for each edge in this entity
413  for( int i = 0; i < num_edges; i++ )
414  {
415  // a node was already assigned
416  if( element[i + num_vertices] != 0 ) continue;
417 
418  tmp_edge_conn[0] = element[entity_edges.conn[i][0]];
419  tmp_edge_conn[1] = element[entity_edges.conn[i][1]];
420 
421  EntityHandle already_made_node = center_node_exist( tmp_edge_conn[0], tmp_edge_conn[1], adjacent_entities );
422 
423  if( already_made_node )
424  {
425  element[i + num_vertices] = already_made_node;
426  }
427  // create a node
428  else
429  {
430  EntitySequence* tmp_sequence = NULL;
431  double sum_coords[3] = { 0, 0, 0 };
432  seq_manager->find( tmp_edge_conn[0], tmp_sequence );
433  static_cast< VertexSequence* >( tmp_sequence )
434  ->get_coordinates( tmp_edge_conn[0], tmp_coords[0], tmp_coords[1], tmp_coords[2] );
435  sum_coords[0] += tmp_coords[0];
436  sum_coords[1] += tmp_coords[1];
437  sum_coords[2] += tmp_coords[2];
438  seq_manager->find( tmp_edge_conn[1], tmp_sequence );
439  static_cast< VertexSequence* >( tmp_sequence )
440  ->get_coordinates( tmp_edge_conn[1], tmp_coords[0], tmp_coords[1], tmp_coords[2] );
441  sum_coords[0] = ( sum_coords[0] + tmp_coords[0] ) / 2;
442  sum_coords[1] = ( sum_coords[1] + tmp_coords[1] ) / 2;
443  sum_coords[2] = ( sum_coords[2] + tmp_coords[2] ) / 2;
444 
445  mMB->create_vertex( sum_coords, element[i + num_vertices] );
446  }
447 
448  if( mHONodeAddedRemoved ) mHONodeAddedRemoved->node_added( element[i + num_vertices], curr_handle );
449  }
450 
451  curr_handle++;
452  }
453 
454  return MB_SUCCESS;
455 }
456 
457 EntityHandle HigherOrderFactory::center_node_exist( EntityHandle corner1,
458  EntityHandle corner2,
459  std::vector< EntityHandle >& adj_entities )
460 {
461  AEntityFactory* a_fact = mMB->a_entity_factory();
462  std::vector< EntityHandle > adj_corner1( 32 );
463  std::vector< EntityHandle > adj_corner2( 32 );
464 
465  // create needed vertex adjacencies
466  if( !a_fact->vert_elem_adjacencies() ) a_fact->create_vert_elem_adjacencies();
467 
468  // vectors are returned sorted
469 
470  a_fact->get_adjacencies( corner1, adj_corner1 );
471  a_fact->get_adjacencies( corner2, adj_corner2 );
472 
473  // these are the entities adjacent to both nodes
474  adj_entities.clear();
475  std::set_intersection( adj_corner1.begin(), adj_corner1.end(), adj_corner2.begin(), adj_corner2.end(),
476  std::back_inserter< std::vector< EntityHandle > >( adj_entities ) );
477 
478  // iterate of the entities to find a mid node
479  const EntityHandle* conn;
480  int conn_size = 0;
481  for( std::vector< EntityHandle >::iterator iter = adj_entities.begin(); iter != adj_entities.end(); )
482  {
483  EntityType this_type = TYPE_FROM_HANDLE( *iter );
484  if( this_type == MBENTITYSET )
485  {
486  ++iter;
487  continue;
488  }
489  mMB->get_connectivity( *iter, conn, conn_size );
490  // if this entity has mid edge nodes
491  if( CN::HasMidEdgeNodes( this_type, conn_size ) )
492  {
493  // find out at which index the mid node should be at
494  int first_node = std::find( conn, conn + conn_size, corner1 ) - conn;
495  int second_node = std::find( conn, conn + conn_size, corner2 ) - conn;
496  if( first_node == conn_size || second_node == conn_size )
497  assert( "We should always find our nodes no matter what" == NULL );
498  int high_node_index = mNodeMap[this_type][first_node][second_node];
499  if( conn[high_node_index] != 0 ) return conn[high_node_index];
500  ++iter;
501  }
502  else
503  {
504  iter = adj_entities.erase( iter );
505  }
506  }
507 
508  return 0;
509 }
510 
511 EntityHandle HigherOrderFactory::center_node_exist( EntityHandle corners[4], std::vector< EntityHandle >& adj_entities )
512 {
513  AEntityFactory* a_fact = mMB->a_entity_factory();
514  std::vector< EntityHandle > adj_corner[4];
515  int num_nodes = corners[3] == 0 ? 3 : 4;
516  int i = 0;
517 
518  // create needed vertex adjacencies
519  if( !a_fact->vert_elem_adjacencies() ) a_fact->create_vert_elem_adjacencies();
520 
521  // vectors are returned sorted
522  for( i = 0; i < num_nodes; i++ )
523  a_fact->get_adjacencies( corners[i], adj_corner[i] );
524 
525  // these are the entities adjacent to both nodes
526  for( i = 1; i < num_nodes; i++ )
527  {
528  adj_entities.clear();
529  std::set_intersection( adj_corner[i - 1].begin(), adj_corner[i - 1].end(), adj_corner[i].begin(),
530  adj_corner[i].end(), std::back_inserter< std::vector< EntityHandle > >( adj_entities ) );
531  adj_corner[i].swap( adj_entities );
532  }
533  adj_entities.swap( adj_corner[i - 1] );
534 
535  // iterate of the entities to find a mid node
536  const EntityHandle* conn;
537  int conn_size = 0;
538  for( std::vector< EntityHandle >::iterator iter = adj_entities.begin(); iter != adj_entities.end(); )
539  {
540  EntityType this_type = TYPE_FROM_HANDLE( *iter );
541  if( this_type == MBENTITYSET )
542  {
543  ++iter;
544  continue;
545  }
546  const CN::ConnMap& entity_faces = CN::mConnectivityMap[this_type][1];
547  mMB->get_connectivity( *iter, conn, conn_size );
548  int offset = CN::VerticesPerEntity( this_type );
549  if( CN::HasMidEdgeNodes( this_type, conn_size ) ) offset += CN::mConnectivityMap[this_type][0].num_sub_elements;
550 
551  // if this entity has mid face nodes
552  if( CN::HasMidFaceNodes( this_type, conn_size ) )
553  {
554  int k;
555  int indexes[4];
556  for( k = 0; k < num_nodes; k++ )
557  indexes[k] = std::find( conn, conn + conn_size, corners[k] ) - conn;
558 
559  // find out at which index the mid node should be at
560  for( k = 0; k < entity_faces.num_sub_elements; k++ )
561  {
562  if( CN::VerticesPerEntity( entity_faces.target_type[k] ) != num_nodes ) continue;
563 
564  int* pivot = std::find( indexes, indexes + num_nodes, entity_faces.conn[k][0] );
565  if( pivot == indexes + num_nodes ) continue;
566 
567  if( pivot != indexes ) std::rotate( indexes, pivot, indexes + num_nodes );
568 
569  if( std::equal( indexes, indexes + num_nodes, entity_faces.conn[k] ) )
570  {
571  if( conn[k + offset] != 0 ) return conn[k + offset];
572  k = entity_faces.num_sub_elements;
573  }
574  else
575  {
576  int temp = indexes[1];
577  indexes[1] = indexes[num_nodes - 1];
578  indexes[num_nodes - 1] = temp;
579  if( std::equal( indexes, indexes + num_nodes, entity_faces.conn[k] ) )
580  {
581  if( conn[k + offset] != 0 ) return conn[k + offset];
582  k = entity_faces.num_sub_elements;
583  }
584  }
585  }
586  ++iter;
587  }
588  else
589  {
590  iter = adj_entities.erase( iter );
591  }
592  }
593 
594  return 0;
595 }
596 
597 bool HigherOrderFactory::add_center_node( EntityType this_type,
598  EntityHandle* element_conn,
599  int conn_size,
600  EntityHandle corner_node1,
601  EntityHandle corner_node2,
602  EntityHandle center_node )
603 {
604  int first_node = std::find( element_conn, element_conn + conn_size, corner_node1 ) - element_conn;
605  int second_node = std::find( element_conn, element_conn + conn_size, corner_node2 ) - element_conn;
606  if( first_node == conn_size || second_node == conn_size )
607  assert( "We should always find our nodes no matter what" == NULL );
608  int high_node_index = mNodeMap[this_type][first_node][second_node];
609  element_conn[high_node_index] = center_node;
610  return true;
611 }
612 
613 ErrorCode HigherOrderFactory::copy_corner_nodes( ElementSequence* src, ElementSequence* dst )
614 {
615  unsigned num_corners = CN::VerticesPerEntity( src->type() );
616  return copy_nodes( src, dst, num_corners, 0, 0 );
617 }
618 
619 ErrorCode HigherOrderFactory::copy_mid_edge_nodes( ElementSequence* src, ElementSequence* dst )
620 {
621  if( !src->has_mid_edge_nodes() || !dst->has_mid_edge_nodes() ) return MB_FAILURE;
622 
623  unsigned num_corners = CN::VerticesPerEntity( src->type() );
624  unsigned num_edges = ( src->type() == MBEDGE ) ? 1 : CN::NumSubEntities( src->type(), 1 );
625  return copy_nodes( src, dst, num_edges, num_corners, num_corners );
626 }
627 
628 ErrorCode HigherOrderFactory::zero_mid_edge_nodes( ElementSequence* dst )
629 {
630  if( !dst->has_mid_edge_nodes() ) return MB_FAILURE;
631 
632  unsigned num_corners = CN::VerticesPerEntity( dst->type() );
633  unsigned num_edges = ( dst->type() == MBEDGE ) ? 1 : CN::NumSubEntities( dst->type(), 1 );
634  return zero_nodes( dst, num_edges, num_corners );
635 }
636 
637 ErrorCode HigherOrderFactory::copy_mid_face_nodes( ElementSequence* src, ElementSequence* dst )
638 {
639  if( !src->has_mid_face_nodes() || !dst->has_mid_face_nodes() ) return MB_FAILURE;
640 
641  unsigned src_offset = CN::VerticesPerEntity( src->type() );
642  unsigned dst_offset = src_offset;
643  if( src->has_mid_edge_nodes() ) src_offset += CN::NumSubEntities( src->type(), 1 );
644  if( dst->has_mid_edge_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 1 );
645  unsigned num_faces = ( CN::Dimension( src->type() ) == 2 ) ? 1 : CN::NumSubEntities( src->type(), 2 );
646  return copy_nodes( src, dst, num_faces, src_offset, dst_offset );
647 }
648 
649 ErrorCode HigherOrderFactory::zero_mid_face_nodes( ElementSequence* dst )
650 {
651  if( !dst->has_mid_face_nodes() ) return MB_FAILURE;
652 
653  unsigned dst_offset = CN::VerticesPerEntity( dst->type() );
654  if( dst->has_mid_edge_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 1 );
655  unsigned num_faces = ( CN::Dimension( dst->type() ) == 2 ) ? 1 : CN::NumSubEntities( dst->type(), 2 );
656  return zero_nodes( dst, num_faces, dst_offset );
657 }
658 
659 ErrorCode HigherOrderFactory::copy_mid_volume_nodes( ElementSequence* src, ElementSequence* dst )
660 {
661  if( !src->has_mid_volume_nodes() || !dst->has_mid_volume_nodes() ) return MB_FAILURE;
662 
663  unsigned src_offset = CN::VerticesPerEntity( src->type() );
664  unsigned dst_offset = src_offset;
665  if( src->has_mid_edge_nodes() ) src_offset += CN::NumSubEntities( src->type(), 1 );
666  if( dst->has_mid_edge_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 1 );
667  if( src->has_mid_face_nodes() ) src_offset += CN::NumSubEntities( src->type(), 2 );
668  if( dst->has_mid_face_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 2 );
669  return copy_nodes( src, dst, 1, src_offset, dst_offset );
670 }
671 
672 ErrorCode HigherOrderFactory::zero_mid_volume_nodes( ElementSequence* dst )
673 {
674  if( !dst->has_mid_volume_nodes() ) return MB_FAILURE;
675 
676  unsigned dst_offset = CN::VerticesPerEntity( dst->type() );
677  if( dst->has_mid_edge_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 1 );
678  if( dst->has_mid_face_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 2 );
679  return zero_nodes( dst, 1, dst_offset );
680 }
681 
682 ErrorCode HigherOrderFactory::copy_nodes( ElementSequence* src,
683  ElementSequence* dst,
684  unsigned nodes_per_elem,
685  unsigned src_offset,
686  unsigned dst_offset )
687 {
688  if( src->type() != dst->type() ) return MB_FAILURE;
689 
690  unsigned src_stride = src->nodes_per_element();
691  unsigned dst_stride = dst->nodes_per_element();
692  EntityHandle* src_conn = src->get_connectivity_array();
693  EntityHandle* dst_conn = dst->get_connectivity_array();
694  if( !src_conn || !dst_conn ) return MB_FAILURE;
695 
696  if( dst->start_handle() < src->start_handle() || dst->end_handle() > src->end_handle() ) return MB_FAILURE;
697 
698  src_conn += ( dst->start_handle() - src->start_handle() ) * src_stride;
699  EntityID count = dst->size();
700  for( EntityID i = 0; i < count; ++i )
701  {
702  for( unsigned j = 0; j < nodes_per_elem; ++j )
703  dst_conn[j + dst_offset] = src_conn[j + src_offset];
704  src_conn += src_stride;
705  dst_conn += dst_stride;
706  }
707 
708  return MB_SUCCESS;
709 }
710 
711 ErrorCode HigherOrderFactory::zero_nodes( ElementSequence* dst, unsigned nodes_per_elem, unsigned offset )
712 {
713  unsigned dst_stride = dst->nodes_per_element();
714  EntityHandle* dst_conn = dst->get_connectivity_array();
715  if( !dst_conn ) return MB_FAILURE;
716 
717  EntityID count = dst->size();
718  for( EntityID i = 0; i < count; ++i )
719  {
720  std::fill( dst_conn + offset, dst_conn + offset + nodes_per_elem, 0 );
721  dst_conn += dst_stride;
722  }
723 
724  return MB_SUCCESS;
725 }
726 
727 ErrorCode HigherOrderFactory::remove_mid_edge_nodes( ElementSequence* seq,
728  EntityHandle start,
729  EntityHandle end,
730  Tag deletable_nodes )
731 {
732  int count;
733  int offset;
734  if( seq->type() == MBEDGE )
735  {
736  count = 1;
737  offset = 2;
738  }
739  else
740  {
741  count = CN::NumSubEntities( seq->type(), 1 );
742  offset = CN::VerticesPerEntity( seq->type() );
743  }
744 
745  return remove_ho_nodes( seq, start, end, count, offset, deletable_nodes );
746 }
747 
748 ErrorCode HigherOrderFactory::remove_mid_face_nodes( ElementSequence* seq,
749  EntityHandle start,
750  EntityHandle end,
751  Tag deletable_nodes )
752 {
753  int count;
754  if( CN::Dimension( seq->type() ) == 2 )
755  count = 1;
756  else
757  count = CN::NumSubEntities( seq->type(), 2 );
758  int offset = CN::VerticesPerEntity( seq->type() );
759  if( seq->has_mid_edge_nodes() ) offset += CN::NumSubEntities( seq->type(), 1 );
760 
761  return remove_ho_nodes( seq, start, end, count, offset, deletable_nodes );
762 }
763 
764 ErrorCode HigherOrderFactory::remove_mid_volume_nodes( ElementSequence* seq,
765  EntityHandle start,
766  EntityHandle end,
767  Tag deletable_nodes )
768 {
769  int offset = CN::VerticesPerEntity( seq->type() );
770  if( seq->has_mid_edge_nodes() ) offset += CN::NumSubEntities( seq->type(), 1 );
771  if( seq->has_mid_face_nodes() ) offset += CN::NumSubEntities( seq->type(), 2 );
772 
773  return remove_ho_nodes( seq, start, end, 1, offset, deletable_nodes );
774 }
775 
776 // Code mostly copied from old EntitySequence.cpp
777 // (ElementEntitySequence::convert_realloc &
778 // ElementEntitySequence::tag_for_deletion).
779 // Copyright from old EntitySequence.cpp:
780 /**
781  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
782  * storing and accessing finite element mesh data.
783  *
784  * Copyright 2004 Sandia Corporation. Under the terms of Contract
785  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
786  * retains certain rights in this software.
787  *
788  * This library is free software; you can redistribute it and/or
789  * modify it under the terms of the GNU Lesser General Public
790  * License as published by the Free Software Foundation; either
791  * version 2.1 of the License, or (at your option) any later version.
792  *
793  */
794 ErrorCode HigherOrderFactory::remove_ho_nodes( ElementSequence* seq,
795  EntityHandle start,
796  EntityHandle end,
797  int nodes_per_elem,
798  int elem_conn_offset,
799  Tag deletable_nodes )
800 {
801  if( start < seq->start_handle() || end > seq->end_handle() ) return MB_ENTITY_NOT_FOUND;
802  EntityHandle* array = seq->get_connectivity_array();
803  if( !array ) return MB_NOT_IMPLEMENTED;
804 
805  std::set< EntityHandle > nodes_processed;
806  for( EntityHandle i = start; i <= end; ++i )
807  { // for each element
808  for( int j = 0; j < nodes_per_elem; ++j )
809  { // for each HO node to remove
810  const EntityID elem = ( i - seq->start_handle() ); // element index
811  const int conn_idx = j + elem_conn_offset;
812  const EntityID index = elem * seq->nodes_per_element() + conn_idx;
813  if( array[index] && nodes_processed.insert( array[index] ).second )
814  {
815  if( tag_for_deletion( i, conn_idx, seq ) )
816  {
817  unsigned char bit = 0x1;
818  mMB->tag_set_data( deletable_nodes, &( array[index] ), 1, &bit );
819  }
820  }
821  }
822  }
823 
824  return MB_SUCCESS;
825 }
826 
827 bool HigherOrderFactory::tag_for_deletion( EntityHandle parent_handle, int conn_index, ElementSequence* seq )
828 {
829  // get type of this sequence
830  EntityType this_type = seq->type();
831 
832  // get dimension of 'parent' element
833  int this_dimension = mMB->dimension_from_handle( parent_handle );
834 
835  // tells us if higher order node is on
836  int dimension, side_number;
837  CN::HONodeParent( this_type, seq->nodes_per_element(), conn_index, dimension, side_number );
838 
839  // it MUST be a higher-order node
840  bool delete_node = false;
841 
842  assert( dimension != -1 );
843  assert( side_number != -1 );
844 
845  // could be a mid-volume/face/edge node on a hex/face/edge respectively
846  // if so...delete it bc/ no one else owns it too
847  std::vector< EntityHandle > connectivity;
848  if( dimension == this_dimension && side_number == 0 )
849  delete_node = true;
850  else // the node could also be on a lower order entity of 'tmp_entity'
851  {
852  // get 'side' of 'parent_handle' that node is on
853  EntityHandle target_entity = 0;
854  mMB->side_element( parent_handle, dimension, side_number, target_entity );
855 
856  if( target_entity )
857  {
858  AEntityFactory* a_fact = mMB->a_entity_factory();
859  EntityHandle low_meshset;
860  int dum;
861  low_meshset = CREATE_HANDLE( MBENTITYSET, 0, dum );
862 
863  // just get corner nodes of target_entity
864  connectivity.clear();
865  ErrorCode rval;
866  rval = mMB->get_connectivity( &( target_entity ), 1, connectivity, true );MB_CHK_ERR( rval );
867 
868  // for each node, get all common adjacencies of nodes in 'parent_handle'
869  std::vector< EntityHandle > adj_list_1, adj_list_2, adj_entities;
870  a_fact->get_adjacencies( connectivity[0], adj_list_1 );
871 
872  // remove meshsets
873  adj_list_1.erase(
874  std::remove_if( adj_list_1.begin(), adj_list_1.end(),
875  std::bind( std::greater< EntityHandle >(), std::placeholders::_1, low_meshset ) ),
876  adj_list_1.end() );
877  // std::bind2nd(std::greater<EntityHandle>(),low_meshset)), adj_list_1.end());
878  // https://stackoverflow.com/questions/32739018/a-replacement-for-stdbind2nd
879 
880  size_t i;
881  for( i = 1; i < connectivity.size(); i++ )
882  {
883  adj_list_2.clear();
884  a_fact->get_adjacencies( connectivity[i], adj_list_2 );
885 
886  // remove meshsets
887  adj_list_2.erase(
888  std::remove_if( adj_list_2.begin(), adj_list_2.end(),
889  std::bind( std::greater< EntityHandle >(), std::placeholders::_1, low_meshset ) ),
890  adj_list_2.end() );
891  // std::bind2nd(std::greater<EntityHandle>(),low_meshset)), adj_list_2.end());
892  // https://stackoverflow.com/questions/32739018/a-replacement-for-stdbind2nd
893 
894  // intersect the 2 lists
895  adj_entities.clear();
896  std::set_intersection( adj_list_1.begin(), adj_list_1.end(), adj_list_2.begin(), adj_list_2.end(),
897  std::back_inserter< std::vector< EntityHandle > >( adj_entities ) );
898  adj_list_1.clear();
899  adj_list_1 = adj_entities;
900  }
901 
902  assert( adj_entities.size() ); // has to have at least one adjacency
903 
904  // see if node is in other elements, not in this sequence...if so, delete it
905  for( i = 0; i < adj_entities.size(); i++ )
906  {
907  if( adj_entities[i] >= seq->start_handle() && adj_entities[i] <= seq->end_handle() )
908  {
909  delete_node = false;
910  break;
911  }
912  else
913  delete_node = true;
914  }
915  }
916  else // there is no lower order entity that also contains node
917  delete_node = true;
918  }
919 
920  return delete_node;
921 }
922 
923 } // namespace moab