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