MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /** 00002 * MOAB, a Mesh-Oriented datABase, is a software component for creating, 00003 * storing and accessing finite element mesh data. 00004 * 00005 * Copyright 2004 Sandia Corporation. Under the terms of Contract 00006 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 00007 * retains certain rights in this software. 00008 * 00009 * This library is free software; you can redistribute it and/or 00010 * modify it under the terms of the GNU Lesser General Public 00011 * License as published by the Free Software Foundation; either 00012 * version 2.1 of the License, or (at your option) any later version. 00013 * 00014 */ 00015 00016 #ifdef WIN32 00017 #ifdef _DEBUG 00018 // turn off warnings that say they debugging identifier has been truncated 00019 // this warning comes up when using some STL containers 00020 #pragma warning( disable : 4786 ) 00021 #endif 00022 #endif 00023 00024 #include "moab/HigherOrderFactory.hpp" 00025 #include "SequenceManager.hpp" 00026 #include "UnstructuredElemSeq.hpp" 00027 #include "VertexSequence.hpp" 00028 #include "AEntityFactory.hpp" 00029 #include "moab/Core.hpp" 00030 #include "moab/CN.hpp" 00031 #include <assert.h> 00032 #include <algorithm> 00033 00034 namespace moab 00035 { 00036 00037 using namespace std; 00038 00039 HigherOrderFactory::HigherOrderFactory( Core* MB, Interface::HONodeAddedRemoved* function_object ) 00040 : mMB( MB ), mHONodeAddedRemoved( function_object ) 00041 { 00042 initialize_map(); 00043 } 00044 HigherOrderFactory::~HigherOrderFactory() {} 00045 00046 // bool HigherOrderFactory::mMapInitialized = false; 00047 00048 void HigherOrderFactory::initialize_map() 00049 { 00050 // if(mMapInitialized) 00051 // return; 00052 00053 for( EntityType i = MBVERTEX; i < MBMAXTYPE; i++ ) 00054 { 00055 const CN::ConnMap& canon_map = CN::mConnectivityMap[i][0]; 00056 unsigned char( &this_map )[8][8] = mNodeMap[i]; 00057 int num_node = CN::VerticesPerEntity( i ); 00058 for( int j = 0; j < canon_map.num_sub_elements; j++ ) 00059 { 00060 unsigned char x = canon_map.conn[j][0]; 00061 unsigned char y = canon_map.conn[j][1]; 00062 this_map[x][y] = num_node; 00063 this_map[y][x] = num_node; 00064 num_node++; 00065 } 00066 } 00067 00068 // mMapInitialized = true; 00069 } 00070 00071 ErrorCode HigherOrderFactory::convert( const EntityHandle meshset, const bool mid_edge_nodes, const bool mid_face_nodes, 00072 const bool mid_volume_nodes ) 00073 { 00074 Range entities; 00075 mMB->get_entities_by_handle( meshset, entities, true ); 00076 return convert( entities, mid_edge_nodes, mid_face_nodes, mid_volume_nodes ); 00077 } 00078 00079 ErrorCode HigherOrderFactory::convert( const Range& entities, const bool mid_edge_nodes, const bool mid_face_nodes, 00080 const bool mid_volume_nodes ) 00081 00082 { 00083 00084 // TODO -- add some more code to prevent from splitting of entity sequences when we don't need 00085 // to. Say we have all hex8's in our mesh and 3 falses are passed in. In the end, no conversion 00086 // will happen, but the sequences could still be split up. 00087 00088 // find out what entity sequences we need to convert 00089 // and where, if necessary, to split them 00090 00091 SequenceManager* seq_manager = mMB->sequence_manager(); 00092 Range::const_pair_iterator p_iter; 00093 for( p_iter = entities.const_pair_begin(); p_iter != entities.const_pair_end(); ++p_iter ) 00094 { 00095 00096 EntityHandle h = p_iter->first; 00097 while( h <= p_iter->second ) 00098 { 00099 00100 EntitySequence* seq; 00101 ErrorCode rval = seq_manager->find( h, seq ); 00102 if( MB_SUCCESS != rval ) return rval; 00103 00104 if( seq->type() == MBVERTEX || seq->type() >= MBENTITYSET ) return MB_TYPE_OUT_OF_RANGE; 00105 00106 // make sequence is not structured mesh 00107 ElementSequence* elemseq = static_cast< ElementSequence* >( seq ); 00108 if( NULL == elemseq->get_connectivity_array() ) return MB_NOT_IMPLEMENTED; 00109 00110 EntityHandle last = p_iter->second; 00111 if( last > seq->end_handle() ) last = seq->end_handle(); 00112 00113 rval = convert_sequence( elemseq, h, last, mid_edge_nodes, mid_face_nodes, mid_volume_nodes ); 00114 if( MB_SUCCESS != rval ) return rval; 00115 00116 h = last + 1; 00117 } 00118 } 00119 00120 return MB_SUCCESS; 00121 } 00122 00123 ErrorCode HigherOrderFactory::convert_sequence( ElementSequence* seq, EntityHandle start, EntityHandle end, 00124 bool mid_edge_nodes, bool mid_face_nodes, bool mid_volume_nodes ) 00125 { 00126 00127 ErrorCode status = MB_SUCCESS; 00128 00129 // lets make sure parameters are ok before we continue 00130 switch( seq->type() ) 00131 { 00132 default: 00133 return MB_TYPE_OUT_OF_RANGE; 00134 case MBEDGE: 00135 mid_face_nodes = false; 00136 case MBTRI: 00137 case MBQUAD: 00138 mid_volume_nodes = false; 00139 case MBTET: 00140 case MBHEX: 00141 case MBPRISM: 00142 case MBPYRAMID: 00143 case MBKNIFE: 00144 break; 00145 } 00146 00147 // calculate number of nodes in target configuration 00148 unsigned nodes_per_elem = CN::VerticesPerEntity( seq->type() ); 00149 if( mid_edge_nodes ) nodes_per_elem += ( seq->type() == MBEDGE ) ? 1 : CN::NumSubEntities( seq->type(), 1 ); 00150 if( mid_face_nodes ) 00151 nodes_per_elem += ( CN::Dimension( seq->type() ) == 2 ) ? 1 : CN::NumSubEntities( seq->type(), 2 ); 00152 if( mid_volume_nodes ) nodes_per_elem += 1; 00153 00154 if( nodes_per_elem == seq->nodes_per_element() ) return MB_SUCCESS; 00155 00156 Tag deletable_nodes; 00157 status = mMB->tag_get_handle( 0, 1, MB_TYPE_BIT, deletable_nodes, MB_TAG_CREAT | MB_TAG_BIT ); 00158 if( MB_SUCCESS != status ) return status; 00159 00160 UnstructuredElemSeq* new_seq = new UnstructuredElemSeq( start, end - start + 1, nodes_per_elem, end - start + 1 ); 00161 00162 copy_corner_nodes( seq, new_seq ); 00163 00164 if( seq->has_mid_edge_nodes() && mid_edge_nodes ) 00165 status = copy_mid_edge_nodes( seq, new_seq ); 00166 else if( seq->has_mid_edge_nodes() && !mid_edge_nodes ) 00167 status = remove_mid_edge_nodes( seq, start, end, deletable_nodes ); 00168 else if( !seq->has_mid_edge_nodes() && mid_edge_nodes ) 00169 status = zero_mid_edge_nodes( new_seq ); 00170 if( MB_SUCCESS != status ) return status; 00171 00172 if( seq->has_mid_face_nodes() && mid_face_nodes ) 00173 status = copy_mid_face_nodes( seq, new_seq ); 00174 else if( seq->has_mid_face_nodes() && !mid_face_nodes ) 00175 status = remove_mid_face_nodes( seq, start, end, deletable_nodes ); 00176 else if( !seq->has_mid_face_nodes() && mid_face_nodes ) 00177 status = zero_mid_face_nodes( new_seq ); 00178 if( MB_SUCCESS != status ) 00179 { 00180 mMB->tag_delete( deletable_nodes ); 00181 return status; 00182 } 00183 00184 if( seq->has_mid_volume_nodes() && mid_volume_nodes ) 00185 status = copy_mid_volume_nodes( seq, new_seq ); 00186 else if( seq->has_mid_volume_nodes() && !mid_volume_nodes ) 00187 status = remove_mid_volume_nodes( seq, start, end, deletable_nodes ); 00188 else if( !seq->has_mid_volume_nodes() && mid_volume_nodes ) 00189 status = zero_mid_volume_nodes( new_seq ); 00190 if( MB_SUCCESS != status ) 00191 { 00192 mMB->tag_delete( deletable_nodes ); 00193 return status; 00194 } 00195 00196 // gather nodes that were marked 00197 Range nodes; 00198 mMB->get_entities_by_type_and_tag( 0, MBVERTEX, &deletable_nodes, NULL, 1, nodes ); 00199 00200 // EntityHandle low_meshset; 00201 // int dum; 00202 // low_meshset = CREATE_HANDLE(MBENTITYSET, 0, dum); 00203 00204 for( Range::iterator iter = nodes.begin(); iter != nodes.end(); ++iter ) 00205 { 00206 unsigned char marked = 0; 00207 mMB->tag_get_data( deletable_nodes, &( *iter ), 1, &marked ); 00208 if( marked ) 00209 { 00210 // we can delete it 00211 if( mHONodeAddedRemoved ) mHONodeAddedRemoved->node_removed( *iter ); 00212 mMB->delete_entities( &( *iter ), 1 ); 00213 } 00214 } 00215 00216 const bool create_midedge = !seq->has_mid_edge_nodes() && mid_edge_nodes; 00217 const bool create_midface = !seq->has_mid_face_nodes() && mid_face_nodes; 00218 const bool create_midvolm = !seq->has_mid_volume_nodes() && mid_volume_nodes; 00219 00220 mMB->tag_delete( deletable_nodes ); 00221 00222 status = mMB->sequence_manager()->replace_subsequence( new_seq ); 00223 if( MB_SUCCESS != status ) 00224 { 00225 SequenceData* data = new_seq->data(); 00226 delete new_seq; 00227 delete data; 00228 return status; 00229 } 00230 00231 if( create_midedge ) 00232 { 00233 status = add_mid_edge_nodes( new_seq ); 00234 if( MB_SUCCESS != status ) return status; 00235 } 00236 if( create_midface ) 00237 { 00238 status = add_mid_face_nodes( new_seq ); 00239 if( MB_SUCCESS != status ) return status; 00240 } 00241 if( create_midvolm ) 00242 { 00243 status = add_mid_volume_nodes( new_seq ); 00244 if( MB_SUCCESS != status ) return status; 00245 } 00246 00247 return status; 00248 } 00249 00250 ErrorCode HigherOrderFactory::add_mid_volume_nodes( ElementSequence* seq ) 00251 { 00252 EntityType this_type = seq->type(); 00253 SequenceManager* seq_manager = mMB->sequence_manager(); 00254 00255 // find out where in the connectivity list to add these new mid volume nodes 00256 int edge_factor = seq->has_mid_edge_nodes() ? 1 : 0; 00257 int face_factor = seq->has_mid_face_nodes() ? 1 : 0; 00258 // offset by number of higher order nodes on edges if they exist 00259 int num_corner_nodes = CN::VerticesPerEntity( this_type ); 00260 int new_node_index = num_corner_nodes; 00261 new_node_index += edge_factor * CN::mConnectivityMap[this_type][0].num_sub_elements; 00262 new_node_index += face_factor * CN::mConnectivityMap[this_type][1].num_sub_elements; 00263 00264 EntityHandle* element = seq->get_connectivity_array(); 00265 EntityHandle curr_handle = seq->start_handle(); 00266 int nodes_per_element = seq->nodes_per_element(); 00267 EntityHandle* end_element = element + nodes_per_element * ( seq->size() ); 00268 00269 // iterate over the elements 00270 for( ; element < end_element; element += nodes_per_element ) 00271 { 00272 // find the centroid of this element 00273 double tmp_coords[3], sum_coords[3] = { 0, 0, 0 }; 00274 EntitySequence* eseq = NULL; 00275 for( int i = 0; i < num_corner_nodes; i++ ) 00276 { 00277 seq_manager->find( element[i], eseq ); 00278 static_cast< VertexSequence* >( eseq )->get_coordinates( element[i], tmp_coords[0], tmp_coords[1], 00279 tmp_coords[2] ); 00280 sum_coords[0] += tmp_coords[0]; 00281 sum_coords[1] += tmp_coords[1]; 00282 sum_coords[2] += tmp_coords[2]; 00283 } 00284 sum_coords[0] /= num_corner_nodes; 00285 sum_coords[1] /= num_corner_nodes; 00286 sum_coords[2] /= num_corner_nodes; 00287 00288 // create a new vertex at the centroid 00289 mMB->create_vertex( sum_coords, element[new_node_index] ); 00290 00291 if( mHONodeAddedRemoved ) mHONodeAddedRemoved->node_added( element[new_node_index], curr_handle ); 00292 00293 curr_handle++; 00294 } 00295 00296 return MB_SUCCESS; 00297 } 00298 00299 ErrorCode HigherOrderFactory::add_mid_face_nodes( ElementSequence* seq ) 00300 { 00301 EntityType this_type = seq->type(); 00302 SequenceManager* seq_manager = mMB->sequence_manager(); 00303 int num_vertices = CN::VerticesPerEntity( this_type ); 00304 int num_edges = CN::mConnectivityMap[this_type][0].num_sub_elements; 00305 num_edges = seq->has_mid_edge_nodes() ? num_edges : 0; 00306 int num_faces = CN::mConnectivityMap[this_type][1].num_sub_elements; 00307 00308 const CN::ConnMap& entity_faces = CN::mConnectivityMap[this_type][1]; 00309 00310 EntityHandle* element = seq->get_connectivity_array(); 00311 EntityHandle curr_handle = seq->start_handle(); 00312 int nodes_per_element = seq->nodes_per_element(); 00313 EntityHandle* end_element = element + nodes_per_element * ( seq->size() ); 00314 00315 EntityHandle tmp_face_conn[4]; // max face nodes = 4 00316 std::vector< EntityHandle > adjacent_entities( 4 ); 00317 00318 double tmp_coords[3]; 00319 00320 // iterate over the elements 00321 for( ; element < end_element; element += nodes_per_element ) 00322 { 00323 // for each edge in this entity 00324 for( int i = 0; i < num_faces; i++ ) 00325 { 00326 // a node was already assigned 00327 if( element[i + num_edges + num_vertices] != 0 ) continue; 00328 00329 tmp_face_conn[0] = element[entity_faces.conn[i][0]]; 00330 tmp_face_conn[1] = element[entity_faces.conn[i][1]]; 00331 tmp_face_conn[2] = element[entity_faces.conn[i][2]]; 00332 if( entity_faces.num_corners_per_sub_element[i] == 4 ) 00333 tmp_face_conn[3] = element[entity_faces.conn[i][3]]; 00334 else 00335 tmp_face_conn[3] = 0; 00336 00337 EntityHandle already_made_node = center_node_exist( tmp_face_conn, adjacent_entities ); 00338 00339 if( already_made_node ) { element[i + num_edges + num_vertices] = already_made_node; } 00340 // create a node 00341 else 00342 { 00343 EntitySequence* tmp_sequence = NULL; 00344 double sum_coords[3] = { 0, 0, 0 }; 00345 int max_nodes = entity_faces.num_corners_per_sub_element[i]; 00346 for( int k = 0; k < max_nodes; k++ ) 00347 { 00348 seq_manager->find( tmp_face_conn[k], tmp_sequence ); 00349 static_cast< VertexSequence* >( tmp_sequence ) 00350 ->get_coordinates( tmp_face_conn[k], tmp_coords[0], tmp_coords[1], tmp_coords[2] ); 00351 sum_coords[0] += tmp_coords[0]; 00352 sum_coords[1] += tmp_coords[1]; 00353 sum_coords[2] += tmp_coords[2]; 00354 } 00355 00356 sum_coords[0] /= max_nodes; 00357 sum_coords[1] /= max_nodes; 00358 sum_coords[2] /= max_nodes; 00359 00360 mMB->create_vertex( sum_coords, element[i + num_edges + num_vertices] ); 00361 } 00362 00363 if( mHONodeAddedRemoved ) 00364 mHONodeAddedRemoved->node_added( element[i + num_edges + num_vertices], curr_handle ); 00365 } 00366 00367 curr_handle++; 00368 } 00369 00370 return MB_SUCCESS; 00371 } 00372 00373 ErrorCode HigherOrderFactory::add_mid_edge_nodes( ElementSequence* seq ) 00374 { 00375 // for each node, need to see if it was already created. 00376 EntityType this_type = seq->type(); 00377 SequenceManager* seq_manager = mMB->sequence_manager(); 00378 00379 // offset by number of corner nodes 00380 int num_vertices = CN::VerticesPerEntity( this_type ); 00381 int num_edges = CN::mConnectivityMap[this_type][0].num_sub_elements; 00382 00383 const CN::ConnMap& entity_edges = CN::mConnectivityMap[this_type][0]; 00384 00385 EntityHandle* element = seq->get_connectivity_array(); 00386 EntityHandle curr_handle = seq->start_handle(); 00387 int nodes_per_element = seq->nodes_per_element(); 00388 EntityHandle* end_element = element + nodes_per_element * ( seq->size() ); 00389 00390 EntityHandle tmp_edge_conn[2]; 00391 std::vector< EntityHandle > adjacent_entities( 32 ); 00392 00393 double tmp_coords[3]; 00394 00395 // iterate over the elements 00396 for( ; element < end_element; element += nodes_per_element ) 00397 { 00398 // for each edge in this entity 00399 for( int i = 0; i < num_edges; i++ ) 00400 { 00401 // a node was already assigned 00402 if( element[i + num_vertices] != 0 ) continue; 00403 00404 tmp_edge_conn[0] = element[entity_edges.conn[i][0]]; 00405 tmp_edge_conn[1] = element[entity_edges.conn[i][1]]; 00406 00407 EntityHandle already_made_node = center_node_exist( tmp_edge_conn[0], tmp_edge_conn[1], adjacent_entities ); 00408 00409 if( already_made_node ) { element[i + num_vertices] = already_made_node; } 00410 // create a node 00411 else 00412 { 00413 EntitySequence* tmp_sequence = NULL; 00414 double sum_coords[3] = { 0, 0, 0 }; 00415 seq_manager->find( tmp_edge_conn[0], tmp_sequence ); 00416 static_cast< VertexSequence* >( tmp_sequence ) 00417 ->get_coordinates( tmp_edge_conn[0], tmp_coords[0], tmp_coords[1], tmp_coords[2] ); 00418 sum_coords[0] += tmp_coords[0]; 00419 sum_coords[1] += tmp_coords[1]; 00420 sum_coords[2] += tmp_coords[2]; 00421 seq_manager->find( tmp_edge_conn[1], tmp_sequence ); 00422 static_cast< VertexSequence* >( tmp_sequence ) 00423 ->get_coordinates( tmp_edge_conn[1], tmp_coords[0], tmp_coords[1], tmp_coords[2] ); 00424 sum_coords[0] = ( sum_coords[0] + tmp_coords[0] ) / 2; 00425 sum_coords[1] = ( sum_coords[1] + tmp_coords[1] ) / 2; 00426 sum_coords[2] = ( sum_coords[2] + tmp_coords[2] ) / 2; 00427 00428 mMB->create_vertex( sum_coords, element[i + num_vertices] ); 00429 } 00430 00431 if( mHONodeAddedRemoved ) mHONodeAddedRemoved->node_added( element[i + num_vertices], curr_handle ); 00432 } 00433 00434 curr_handle++; 00435 } 00436 00437 return MB_SUCCESS; 00438 } 00439 00440 EntityHandle HigherOrderFactory::center_node_exist( EntityHandle corner1, EntityHandle corner2, 00441 std::vector< EntityHandle >& adj_entities ) 00442 { 00443 AEntityFactory* a_fact = mMB->a_entity_factory(); 00444 std::vector< EntityHandle > adj_corner1( 32 ); 00445 std::vector< EntityHandle > adj_corner2( 32 ); 00446 00447 // create needed vertex adjacencies 00448 if( !a_fact->vert_elem_adjacencies() ) a_fact->create_vert_elem_adjacencies(); 00449 00450 // vectors are returned sorted 00451 00452 a_fact->get_adjacencies( corner1, adj_corner1 ); 00453 a_fact->get_adjacencies( corner2, adj_corner2 ); 00454 00455 // these are the entities adjacent to both nodes 00456 adj_entities.clear(); 00457 std::set_intersection( adj_corner1.begin(), adj_corner1.end(), adj_corner2.begin(), adj_corner2.end(), 00458 std::back_inserter< std::vector< EntityHandle > >( adj_entities ) ); 00459 00460 // iterate of the entities to find a mid node 00461 const EntityHandle* conn; 00462 int conn_size = 0; 00463 for( std::vector< EntityHandle >::iterator iter = adj_entities.begin(); iter != adj_entities.end(); ) 00464 { 00465 EntityType this_type = TYPE_FROM_HANDLE( *iter ); 00466 if( this_type == MBENTITYSET ) 00467 { 00468 ++iter; 00469 continue; 00470 } 00471 mMB->get_connectivity( *iter, conn, conn_size ); 00472 // if this entity has mid edge nodes 00473 if( CN::HasMidEdgeNodes( this_type, conn_size ) ) 00474 { 00475 // find out at which index the mid node should be at 00476 int first_node = std::find( conn, conn + conn_size, corner1 ) - conn; 00477 int second_node = std::find( conn, conn + conn_size, corner2 ) - conn; 00478 if( first_node == conn_size || second_node == conn_size ) 00479 assert( "We should always find our nodes no matter what" == NULL ); 00480 int high_node_index = mNodeMap[this_type][first_node][second_node]; 00481 if( conn[high_node_index] != 0 ) return conn[high_node_index]; 00482 ++iter; 00483 } 00484 else 00485 { 00486 iter = adj_entities.erase( iter ); 00487 } 00488 } 00489 00490 return 0; 00491 } 00492 00493 EntityHandle HigherOrderFactory::center_node_exist( EntityHandle corners[4], std::vector< EntityHandle >& adj_entities ) 00494 { 00495 AEntityFactory* a_fact = mMB->a_entity_factory(); 00496 std::vector< EntityHandle > adj_corner[4]; 00497 int num_nodes = corners[3] == 0 ? 3 : 4; 00498 int i = 0; 00499 00500 // create needed vertex adjacencies 00501 if( !a_fact->vert_elem_adjacencies() ) a_fact->create_vert_elem_adjacencies(); 00502 00503 // vectors are returned sorted 00504 for( i = 0; i < num_nodes; i++ ) 00505 a_fact->get_adjacencies( corners[i], adj_corner[i] ); 00506 00507 // these are the entities adjacent to both nodes 00508 for( i = 1; i < num_nodes; i++ ) 00509 { 00510 adj_entities.clear(); 00511 std::set_intersection( adj_corner[i - 1].begin(), adj_corner[i - 1].end(), adj_corner[i].begin(), 00512 adj_corner[i].end(), std::back_inserter< std::vector< EntityHandle > >( adj_entities ) ); 00513 adj_corner[i].swap( adj_entities ); 00514 } 00515 adj_entities.swap( adj_corner[i - 1] ); 00516 00517 // iterate of the entities to find a mid node 00518 const EntityHandle* conn; 00519 int conn_size = 0; 00520 for( std::vector< EntityHandle >::iterator iter = adj_entities.begin(); iter != adj_entities.end(); ) 00521 { 00522 EntityType this_type = TYPE_FROM_HANDLE( *iter ); 00523 if( this_type == MBENTITYSET ) 00524 { 00525 ++iter; 00526 continue; 00527 } 00528 const CN::ConnMap& entity_faces = CN::mConnectivityMap[this_type][1]; 00529 mMB->get_connectivity( *iter, conn, conn_size ); 00530 int offset = CN::VerticesPerEntity( this_type ); 00531 if( CN::HasMidEdgeNodes( this_type, conn_size ) ) offset += CN::mConnectivityMap[this_type][0].num_sub_elements; 00532 00533 // if this entity has mid face nodes 00534 if( CN::HasMidFaceNodes( this_type, conn_size ) ) 00535 { 00536 int k; 00537 int indexes[4]; 00538 for( k = 0; k < num_nodes; k++ ) 00539 indexes[k] = std::find( conn, conn + conn_size, corners[k] ) - conn; 00540 00541 // find out at which index the mid node should be at 00542 for( k = 0; k < entity_faces.num_sub_elements; k++ ) 00543 { 00544 if( CN::VerticesPerEntity( entity_faces.target_type[k] ) != num_nodes ) continue; 00545 00546 int* pivot = std::find( indexes, indexes + num_nodes, entity_faces.conn[k][0] ); 00547 if( pivot == indexes + num_nodes ) continue; 00548 00549 if( pivot != indexes ) std::rotate( indexes, pivot, indexes + num_nodes ); 00550 00551 if( std::equal( indexes, indexes + num_nodes, entity_faces.conn[k] ) ) 00552 { 00553 if( conn[k + offset] != 0 ) return conn[k + offset]; 00554 k = entity_faces.num_sub_elements; 00555 } 00556 else 00557 { 00558 int temp = indexes[1]; 00559 indexes[1] = indexes[num_nodes - 1]; 00560 indexes[num_nodes - 1] = temp; 00561 if( std::equal( indexes, indexes + num_nodes, entity_faces.conn[k] ) ) 00562 { 00563 if( conn[k + offset] != 0 ) return conn[k + offset]; 00564 k = entity_faces.num_sub_elements; 00565 } 00566 } 00567 } 00568 ++iter; 00569 } 00570 else 00571 { 00572 iter = adj_entities.erase( iter ); 00573 } 00574 } 00575 00576 return 0; 00577 } 00578 00579 bool HigherOrderFactory::add_center_node( EntityType this_type, EntityHandle* element_conn, int conn_size, 00580 EntityHandle corner_node1, EntityHandle corner_node2, 00581 EntityHandle center_node ) 00582 { 00583 int first_node = std::find( element_conn, element_conn + conn_size, corner_node1 ) - element_conn; 00584 int second_node = std::find( element_conn, element_conn + conn_size, corner_node2 ) - element_conn; 00585 if( first_node == conn_size || second_node == conn_size ) 00586 assert( "We should always find our nodes no matter what" == NULL ); 00587 int high_node_index = mNodeMap[this_type][first_node][second_node]; 00588 element_conn[high_node_index] = center_node; 00589 return true; 00590 } 00591 00592 ErrorCode HigherOrderFactory::copy_corner_nodes( ElementSequence* src, ElementSequence* dst ) 00593 { 00594 unsigned num_corners = CN::VerticesPerEntity( src->type() ); 00595 return copy_nodes( src, dst, num_corners, 0, 0 ); 00596 } 00597 00598 ErrorCode HigherOrderFactory::copy_mid_edge_nodes( ElementSequence* src, ElementSequence* dst ) 00599 { 00600 if( !src->has_mid_edge_nodes() || !dst->has_mid_edge_nodes() ) return MB_FAILURE; 00601 00602 unsigned num_corners = CN::VerticesPerEntity( src->type() ); 00603 unsigned num_edges = ( src->type() == MBEDGE ) ? 1 : CN::NumSubEntities( src->type(), 1 ); 00604 return copy_nodes( src, dst, num_edges, num_corners, num_corners ); 00605 } 00606 00607 ErrorCode HigherOrderFactory::zero_mid_edge_nodes( ElementSequence* dst ) 00608 { 00609 if( !dst->has_mid_edge_nodes() ) return MB_FAILURE; 00610 00611 unsigned num_corners = CN::VerticesPerEntity( dst->type() ); 00612 unsigned num_edges = ( dst->type() == MBEDGE ) ? 1 : CN::NumSubEntities( dst->type(), 1 ); 00613 return zero_nodes( dst, num_edges, num_corners ); 00614 } 00615 00616 ErrorCode HigherOrderFactory::copy_mid_face_nodes( ElementSequence* src, ElementSequence* dst ) 00617 { 00618 if( !src->has_mid_face_nodes() || !dst->has_mid_face_nodes() ) return MB_FAILURE; 00619 00620 unsigned src_offset = CN::VerticesPerEntity( src->type() ); 00621 unsigned dst_offset = src_offset; 00622 if( src->has_mid_edge_nodes() ) src_offset += CN::NumSubEntities( src->type(), 1 ); 00623 if( dst->has_mid_edge_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 1 ); 00624 unsigned num_faces = ( CN::Dimension( src->type() ) == 2 ) ? 1 : CN::NumSubEntities( src->type(), 2 ); 00625 return copy_nodes( src, dst, num_faces, src_offset, dst_offset ); 00626 } 00627 00628 ErrorCode HigherOrderFactory::zero_mid_face_nodes( ElementSequence* dst ) 00629 { 00630 if( !dst->has_mid_face_nodes() ) return MB_FAILURE; 00631 00632 unsigned dst_offset = CN::VerticesPerEntity( dst->type() ); 00633 if( dst->has_mid_edge_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 1 ); 00634 unsigned num_faces = ( CN::Dimension( dst->type() ) == 2 ) ? 1 : CN::NumSubEntities( dst->type(), 2 ); 00635 return zero_nodes( dst, num_faces, dst_offset ); 00636 } 00637 00638 ErrorCode HigherOrderFactory::copy_mid_volume_nodes( ElementSequence* src, ElementSequence* dst ) 00639 { 00640 if( !src->has_mid_volume_nodes() || !dst->has_mid_volume_nodes() ) return MB_FAILURE; 00641 00642 unsigned src_offset = CN::VerticesPerEntity( src->type() ); 00643 unsigned dst_offset = src_offset; 00644 if( src->has_mid_edge_nodes() ) src_offset += CN::NumSubEntities( src->type(), 1 ); 00645 if( dst->has_mid_edge_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 1 ); 00646 if( src->has_mid_face_nodes() ) src_offset += CN::NumSubEntities( src->type(), 2 ); 00647 if( dst->has_mid_face_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 2 ); 00648 return copy_nodes( src, dst, 1, src_offset, dst_offset ); 00649 } 00650 00651 ErrorCode HigherOrderFactory::zero_mid_volume_nodes( ElementSequence* dst ) 00652 { 00653 if( !dst->has_mid_volume_nodes() ) return MB_FAILURE; 00654 00655 unsigned dst_offset = CN::VerticesPerEntity( dst->type() ); 00656 if( dst->has_mid_edge_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 1 ); 00657 if( dst->has_mid_face_nodes() ) dst_offset += CN::NumSubEntities( dst->type(), 2 ); 00658 return zero_nodes( dst, 1, dst_offset ); 00659 } 00660 00661 ErrorCode HigherOrderFactory::copy_nodes( ElementSequence* src, ElementSequence* dst, unsigned nodes_per_elem, 00662 unsigned src_offset, unsigned dst_offset ) 00663 { 00664 if( src->type() != dst->type() ) return MB_FAILURE; 00665 00666 unsigned src_stride = src->nodes_per_element(); 00667 unsigned dst_stride = dst->nodes_per_element(); 00668 EntityHandle* src_conn = src->get_connectivity_array(); 00669 EntityHandle* dst_conn = dst->get_connectivity_array(); 00670 if( !src_conn || !dst_conn ) return MB_FAILURE; 00671 00672 if( dst->start_handle() < src->start_handle() || dst->end_handle() > src->end_handle() ) return MB_FAILURE; 00673 00674 src_conn += ( dst->start_handle() - src->start_handle() ) * src_stride; 00675 EntityID count = dst->size(); 00676 for( EntityID i = 0; i < count; ++i ) 00677 { 00678 for( unsigned j = 0; j < nodes_per_elem; ++j ) 00679 dst_conn[j + dst_offset] = src_conn[j + src_offset]; 00680 src_conn += src_stride; 00681 dst_conn += dst_stride; 00682 } 00683 00684 return MB_SUCCESS; 00685 } 00686 00687 ErrorCode HigherOrderFactory::zero_nodes( ElementSequence* dst, unsigned nodes_per_elem, unsigned offset ) 00688 { 00689 unsigned dst_stride = dst->nodes_per_element(); 00690 EntityHandle* dst_conn = dst->get_connectivity_array(); 00691 if( !dst_conn ) return MB_FAILURE; 00692 00693 EntityID count = dst->size(); 00694 for( EntityID i = 0; i < count; ++i ) 00695 { 00696 std::fill( dst_conn + offset, dst_conn + offset + nodes_per_elem, 0 ); 00697 dst_conn += dst_stride; 00698 } 00699 00700 return MB_SUCCESS; 00701 } 00702 00703 ErrorCode HigherOrderFactory::remove_mid_edge_nodes( ElementSequence* seq, EntityHandle start, EntityHandle end, 00704 Tag deletable_nodes ) 00705 { 00706 int count; 00707 int offset; 00708 if( seq->type() == MBEDGE ) 00709 { 00710 count = 1; 00711 offset = 2; 00712 } 00713 else 00714 { 00715 count = CN::NumSubEntities( seq->type(), 1 ); 00716 offset = CN::VerticesPerEntity( seq->type() ); 00717 } 00718 00719 return remove_ho_nodes( seq, start, end, count, offset, deletable_nodes ); 00720 } 00721 00722 ErrorCode HigherOrderFactory::remove_mid_face_nodes( ElementSequence* seq, EntityHandle start, EntityHandle end, 00723 Tag deletable_nodes ) 00724 { 00725 int count; 00726 if( CN::Dimension( seq->type() ) == 2 ) 00727 count = 1; 00728 else 00729 count = CN::NumSubEntities( seq->type(), 2 ); 00730 int offset = CN::VerticesPerEntity( seq->type() ); 00731 if( seq->has_mid_edge_nodes() ) offset += CN::NumSubEntities( seq->type(), 1 ); 00732 00733 return remove_ho_nodes( seq, start, end, count, offset, deletable_nodes ); 00734 } 00735 00736 ErrorCode HigherOrderFactory::remove_mid_volume_nodes( ElementSequence* seq, EntityHandle start, EntityHandle end, 00737 Tag deletable_nodes ) 00738 { 00739 int offset = CN::VerticesPerEntity( seq->type() ); 00740 if( seq->has_mid_edge_nodes() ) offset += CN::NumSubEntities( seq->type(), 1 ); 00741 if( seq->has_mid_face_nodes() ) offset += CN::NumSubEntities( seq->type(), 2 ); 00742 00743 return remove_ho_nodes( seq, start, end, 1, offset, deletable_nodes ); 00744 } 00745 00746 // Code mostly copied from old EntitySequence.cpp 00747 // (ElementEntitySequence::convert_realloc & 00748 // ElementEntitySequence::tag_for_deletion). 00749 // Copyright from old EntitySequence.cpp: 00750 /** 00751 * MOAB, a Mesh-Oriented datABase, is a software component for creating, 00752 * storing and accessing finite element mesh data. 00753 * 00754 * Copyright 2004 Sandia Corporation. Under the terms of Contract 00755 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 00756 * retains certain rights in this software. 00757 * 00758 * This library is free software; you can redistribute it and/or 00759 * modify it under the terms of the GNU Lesser General Public 00760 * License as published by the Free Software Foundation; either 00761 * version 2.1 of the License, or (at your option) any later version. 00762 * 00763 */ 00764 ErrorCode HigherOrderFactory::remove_ho_nodes( ElementSequence* seq, EntityHandle start, EntityHandle end, 00765 int nodes_per_elem, int elem_conn_offset, Tag deletable_nodes ) 00766 { 00767 if( start < seq->start_handle() || end > seq->end_handle() ) return MB_ENTITY_NOT_FOUND; 00768 EntityHandle* array = seq->get_connectivity_array(); 00769 if( !array ) return MB_NOT_IMPLEMENTED; 00770 00771 std::set< EntityHandle > nodes_processed; 00772 for( EntityHandle i = start; i <= end; ++i ) 00773 { // for each element 00774 for( int j = 0; j < nodes_per_elem; ++j ) 00775 { // for each HO node to remove 00776 const EntityID elem = ( i - seq->start_handle() ); // element index 00777 const int conn_idx = j + elem_conn_offset; 00778 const EntityID index = elem * seq->nodes_per_element() + conn_idx; 00779 if( array[index] && nodes_processed.insert( array[index] ).second ) 00780 { 00781 if( tag_for_deletion( i, conn_idx, seq ) ) 00782 { 00783 unsigned char bit = 0x1; 00784 mMB->tag_set_data( deletable_nodes, &( array[index] ), 1, &bit ); 00785 } 00786 } 00787 } 00788 } 00789 00790 return MB_SUCCESS; 00791 } 00792 00793 bool HigherOrderFactory::tag_for_deletion( EntityHandle parent_handle, int conn_index, ElementSequence* seq ) 00794 { 00795 // get type of this sequence 00796 EntityType this_type = seq->type(); 00797 00798 // get dimension of 'parent' element 00799 int this_dimension = mMB->dimension_from_handle( parent_handle ); 00800 00801 // tells us if higher order node is on 00802 int dimension, side_number; 00803 CN::HONodeParent( this_type, seq->nodes_per_element(), conn_index, dimension, side_number ); 00804 00805 // it MUST be a higher-order node 00806 bool delete_node = false; 00807 00808 assert( dimension != -1 ); 00809 assert( side_number != -1 ); 00810 00811 // could be a mid-volume/face/edge node on a hex/face/edge respectively 00812 // if so...delete it bc/ no one else owns it too 00813 std::vector< EntityHandle > connectivity; 00814 if( dimension == this_dimension && side_number == 0 ) 00815 delete_node = true; 00816 else // the node could also be on a lower order entity of 'tmp_entity' 00817 { 00818 // get 'side' of 'parent_handle' that node is on 00819 EntityHandle target_entity = 0; 00820 mMB->side_element( parent_handle, dimension, side_number, target_entity ); 00821 00822 if( target_entity ) 00823 { 00824 AEntityFactory* a_fact = mMB->a_entity_factory(); 00825 EntityHandle low_meshset; 00826 int dum; 00827 low_meshset = CREATE_HANDLE( MBENTITYSET, 0, dum ); 00828 00829 // just get corner nodes of target_entity 00830 connectivity.clear(); 00831 ErrorCode rval; 00832 rval = mMB->get_connectivity( &( target_entity ), 1, connectivity, true );MB_CHK_ERR( rval ); 00833 00834 // for each node, get all common adjacencies of nodes in 'parent_handle' 00835 std::vector< EntityHandle > adj_list_1, adj_list_2, adj_entities; 00836 a_fact->get_adjacencies( connectivity[0], adj_list_1 ); 00837 00838 // remove meshsets 00839 adj_list_1.erase( 00840 std::remove_if( adj_list_1.begin(), adj_list_1.end(), 00841 std::bind( std::greater< EntityHandle >(), std::placeholders::_1, low_meshset ) ), 00842 adj_list_1.end() ); 00843 // std::bind2nd(std::greater<EntityHandle>(),low_meshset)), adj_list_1.end()); 00844 // https://stackoverflow.com/questions/32739018/a-replacement-for-stdbind2nd 00845 00846 size_t i; 00847 for( i = 1; i < connectivity.size(); i++ ) 00848 { 00849 adj_list_2.clear(); 00850 a_fact->get_adjacencies( connectivity[i], adj_list_2 ); 00851 00852 // remove meshsets 00853 adj_list_2.erase( 00854 std::remove_if( adj_list_2.begin(), adj_list_2.end(), 00855 std::bind( std::greater< EntityHandle >(), std::placeholders::_1, low_meshset ) ), 00856 adj_list_2.end() ); 00857 // std::bind2nd(std::greater<EntityHandle>(),low_meshset)), adj_list_2.end()); 00858 // https://stackoverflow.com/questions/32739018/a-replacement-for-stdbind2nd 00859 00860 // intersect the 2 lists 00861 adj_entities.clear(); 00862 std::set_intersection( adj_list_1.begin(), adj_list_1.end(), adj_list_2.begin(), adj_list_2.end(), 00863 std::back_inserter< std::vector< EntityHandle > >( adj_entities ) ); 00864 adj_list_1.clear(); 00865 adj_list_1 = adj_entities; 00866 } 00867 00868 assert( adj_entities.size() ); // has to have at least one adjacency 00869 00870 // see if node is in other elements, not in this sequence...if so, delete it 00871 for( i = 0; i < adj_entities.size(); i++ ) 00872 { 00873 if( adj_entities[i] >= seq->start_handle() && adj_entities[i] <= seq->end_handle() ) 00874 { 00875 delete_node = false; 00876 break; 00877 } 00878 else 00879 delete_node = true; 00880 } 00881 } 00882 else // there is no lower order entity that also contains node 00883 delete_node = true; 00884 } 00885 00886 return delete_node; 00887 } 00888 00889 } // namespace moab