FaceSets.h

Go to the documentation of this file.

#ifndef __smoab_FaceSets_h
#define __smoab_FaceSets_h

#include "CellSets.h"
#include <set>

namespace smoab
{
//----------------------------------------------------------------------------
class FaceCellSet : public CellSet
{
  public:
    FaceCellSet( int id, smoab::EntityHandle p, const smoab::Range& cells ) : CellSet( p, cells ), ID( id ) {}

    int faceId() const
    {
        return ID;
    }
    void overrideFaceId( int i )
    {
        ID = i;
    }  //USE AT YOUR OWN RISK

  private:
    int ID;
};

//----------------------------------------------------------------------------
typedef std::vector< FaceCellSet > FaceCellSets;

//----------------------------------------------------------------------------
//class that store the regions that a faces is adjacent too
struct FacesAdjRegions
{
    FacesAdjRegions( int f, smoab::EntityHandle r0, smoab::EntityHandle r1 ) : FaceId( f ), Region0( r0 ), Region1( r1 )
    {
        if( r0 > r1 )
        {
            std::swap( this->Region0, this->Region1 );
        }
    }

    FacesAdjRegions( int f ) : FaceId( f ), Region0( -3 ), Region1( -2 ) {}

    bool operator<( const FacesAdjRegions& other ) const
    {
        return ( this->FaceId < other.FaceId );
    }

    smoab::EntityHandle otherId( smoab::EntityHandle other ) const
    {
        if( other == Region0 )
        {
            return Region1;
        }
        return Region0;
    }

    int FaceId;
    smoab::EntityHandle Region0;
    smoab::EntityHandle Region1;
};
//----------------------------------------------------------------------------
smoab::FaceCellSets findFaceSets( smoab::CellSets shells,
                                  smoab::CellSets boundaries,
                                  std::set< smoab::FacesAdjRegions >& faceMaps )
{
    typedef smoab::CellSets::iterator iterator;
    typedef smoab::FaceCellSets::iterator faceIterator;
    typedef std::set< smoab::FacesAdjRegions >::const_iterator FaceAdjIterator;

    //we need to properly label each unique face in shells
    //we do this by intersecting each shell with each other shell
    //to find shell on shell contact, and than we intersect each
    //resulting shell with the boundary conditions
    //the end result of these intersections will be the new modelfaces
    int faceId = 1;
    smoab::FaceCellSets shellFaces;

    //first intersect each shell with each other shell
    std::set< smoab::FacesAdjRegions > shellFaceContacts;
    for( iterator i = shells.begin(); i != shells.end(); ++i )
    {
        //copy the cells so we can add a face that represents
        //all the cells of the region that aren't shared with another region
        int numCells = i->cells().size();  //size() on range is slow, so cache it
        for( iterator j = i + 1; j != shells.end() && numCells > 0; ++j )
        {
            //intersect i and j to make a new face
            smoab::Range intersection = smoab::intersect( i->cells(), j->cells() );
            if( !intersection.empty() )
            {
                //don't want to increment faceId when the intersection is empty
                smoab::FaceCellSet face( faceId++, i->entity(), intersection );
                shellFaces.push_back( face );
                i->erase( intersection );
                j->erase( intersection );
                numCells -= intersection.size();

                //add this to the face map
                smoab::FacesAdjRegions faceInfo( faceId - 1, i->entity(), j->entity() );
                shellFaceContacts.insert( faceInfo );
            }
        }
        //if all the cells for shell i are used, don't add a new
        //empty face
        if( numCells > 0 )
        {
            smoab::FaceCellSet face( faceId++, i->entity(), i->cells() );
            shellFaces.push_back( face );

            //add this to the face map
            smoab::FacesAdjRegions faceInfo( faceId - 1, -1, i->entity() );
            shellFaceContacts.insert( faceInfo );
        }
    }

    //now we have all the faces that match shell on shell contact
    //we know process all the new faces to see if they intersect
    //with any boundary sets. A boundary set can span multiple
    //shells so we want to process it as a second loop

    //store the end before we start adding boundary faces, which
    //we don't need to check agianst other boundaries
    faceId = 1;  //reset the faced id

    //store in a new face set, expanding the current one causes incorrect results
    smoab::FaceCellSets faces;
    for( faceIterator i = shellFaces.begin(); i != shellFaces.end(); ++i )
    {
        //determine from the shell faces if the new face we are creating
        //is bounded by two regions or just one
        smoab::FacesAdjRegions idToSearchFor( i->faceId() );
        FaceAdjIterator adjRegions        = shellFaceContacts.find( idToSearchFor );
        smoab::EntityHandle otherRegionId = adjRegions->otherId( i->entity() );

        int numCells = i->cells().size();  //size() on range is slow, so cache it
        for( iterator j = boundaries.begin(); j != boundaries.end(); ++j )
        {
            smoab::Range intersect = smoab::intersect( i->cells(), j->cells() );
            if( !intersect.empty() )
            {
                //don't want to increment faceId when the intersection is empty
                smoab::FaceCellSet face( faceId++, j->entity(), intersect );
                faces.push_back( face );
                i->erase( intersect );
                numCells -= intersect.size();
                smoab::FacesAdjRegions faceInfo( faceId - 1, i->entity(), otherRegionId );
                faceMaps.insert( faceInfo );
            }
        }
        if( numCells > 0 )
        {
            smoab::FaceCellSet face( faceId++, i->entity(), i->cells() );
            faces.push_back( face );
            smoab::FacesAdjRegions faceInfo( faceId - 1, i->entity(), otherRegionId );
            faceMaps.insert( faceInfo );
        }
    }
    return faces;
}

//----------------------------------------------------------------------------
template < typename T >
std::vector< T > faceIdsPerCell( const smoab::FaceCellSets& faces )
{
    typedef smoab::FaceCellSets::const_iterator iterator;
    typedef std::vector< smoab::EntityHandle >::const_iterator EntityHandleIterator;
    typedef smoab::Range::const_iterator RangeIterator;

    //find all the cells that are in the faceCellSet, and than map
    //the proper face id to that relative position, here comes lower_bounds!
    std::vector< smoab::EntityHandle > searchableCells;
    smoab::Range faceRange = smoab::getAllCells( faces );
    smoab::RangeToVector( faceRange, searchableCells );

    //faceIds will be the resulting array
    std::vector< T > faceIds( searchableCells.size() );

    //construct the start and end iterators for the lower bounds call

    EntityHandleIterator s_begin = searchableCells.begin();
    EntityHandleIterator s_end   = searchableCells.end();

    //search the face cell sets
    for( iterator i = faces.begin(); i != faces.end(); ++i )
    {
        T value                           = static_cast< T >( i->faceId() );
        const smoab::Range& entitiesCells = i->cells();
        for( RangeIterator j = entitiesCells.begin(); j != entitiesCells.end(); ++j )
        {
            EntityHandleIterator result = std::lower_bound( s_begin, s_end, *j );
            std::size_t newId           = std::distance( s_begin, result );
            faceIds[newId]              = value;
        }
    }
    return faceIds;
}

//----------------------------------------------------------------------------
//given a face adjacency, determine the regions spare tag values
template < typename T >
std::pair< T, T > FaceAdjRegionValues( const smoab::FacesAdjRegions& faceAdj,
                                       smoab::Tag* t,
                                       const smoab::Interface& interface )
{
    std::pair< T, T > returnValue;
    const int defaultValue = interface.getDefaultTagVaue< int >( *t );

    /*
   *  IF A REGION IS SET TO -1 WE NEED TO PUSH THAT VALUE DOWN
   *  AS THE MATERIAL, SINCE THE MOAB DEFAULT TAG VALUE WILL
   *  BE CONSIDIERED A REGION, AND WE WANT TO SAY IT BOUNDS THE
   *  VOID REGION
   */
    int tagValue = defaultValue;  //use tagValue to pass in default value
    if( faceAdj.Region0 != -1 )
    {
        tagValue = interface.getTagData( *t, faceAdj.Region0, tagValue );
    }
    else
    {
        tagValue = -1;
    }

    //set the first region tag value into the pair we are returing
    returnValue.first = static_cast< T >( tagValue );

    tagValue = defaultValue;  //use tagValue to pass in default value
    tagValue = interface.getTagData( *t, faceAdj.Region1, tagValue );
    if( faceAdj.Region1 != -1 )
    {
        tagValue = interface.getTagData( *t, faceAdj.Region1, tagValue );
    }
    else
    {
        tagValue = -1;
    }
    returnValue.second = static_cast< T >( tagValue );

    return returnValue;
}

}  // namespace smoab

#endif