MOAB: Mesh Oriented datABase  (version 5.2.1)
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 <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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines