MOAB: Mesh Oriented datABase  (version 5.2.1)
moab::LloydSmoother Class Reference

#include <LloydSmoother.hpp>

+ Collaboration diagram for moab::LloydSmoother:

Public Member Functions

 LloydSmoother (Interface *impl)
 LloydSmoother (Interface *impl, ParallelComm *pc, Range &elems, Tag cds_tag=0, Tag fixed_tag=0, double abs_tol=-1.0, double rel_tol=1.0e-6)
 ~LloydSmoother ()
ErrorCode perform_smooth ()
Interface * mb_impl ()
ParallelCommpcomm ()
void pcomm (ParallelComm *pc)
Rangeelems ()
const Rangeelems () const
Tag fixed_tag ()
void fixed_tag (Tag fixed)
Tag coords_tag ()
void coords_tag (Tag coords)
double abs_tol ()
void abs_tol (double tol)
double rel_tol ()
void rel_tol (double tol)
int num_its ()
void num_its (int num)
int report_its ()
void report_its (int num)

Private Member Functions

ErrorCode initialize ()

Private Attributes

Interface * mbImpl
ParallelCommmyPcomm
Range myElems
Tag coordsTag
Tag fixedTag
double absTol
double relTol
int reportIts
int numIts
bool iCreatedTag

Detailed Description

Examples:
DeformMeshRemap.cpp.

Definition at line 26 of file LloydSmoother.hpp.


Constructor & Destructor Documentation

moab::LloydSmoother::LloydSmoother ( Interface *  impl)
moab::LloydSmoother::LloydSmoother ( Interface *  impl,
ParallelComm pc,
Range elems,
Tag  cds_tag = 0,
Tag  fixed_tag = 0,
double  abs_tol = -1.0,
double  rel_tol = 1.0e-6 
)

Definition at line 17 of file LloydSmoother.cpp.

    : mbImpl( impl ), myPcomm( pc ), myElems( elms ), coordsTag( ctag ), fixedTag( ftag ), absTol( at ), relTol( rt ),
      reportIts( 0 ), numIts( 0 ), iCreatedTag( false )
{
}

Definition at line 23 of file LloydSmoother.cpp.

References fixedTag, iCreatedTag, MB_CHK_SET_ERR_RET, and mbImpl.

{
    if( iCreatedTag && fixedTag )
    {
        ErrorCode rval = mbImpl->tag_delete( fixedTag );MB_CHK_SET_ERR_RET( rval, "Failed to delete the fixed tag" );
    }
}

Member Function Documentation

double moab::LloydSmoother::abs_tol ( ) [inline]

Definition at line 122 of file LloydSmoother.hpp.

References absTol.

    {
        return absTol;
    }
void moab::LloydSmoother::abs_tol ( double  tol) [inline]

Definition at line 129 of file LloydSmoother.hpp.

References absTol.

    {
        absTol = tol;
    }

Definition at line 108 of file LloydSmoother.hpp.

References coordsTag.

    {
        return coordsTag;
    }
void moab::LloydSmoother::coords_tag ( Tag  coords) [inline]

Definition at line 115 of file LloydSmoother.hpp.

References coordsTag.

    {
        coordsTag = coords;
    }

Definition at line 80 of file LloydSmoother.hpp.

References myElems.

    {
        return myElems;
    }
const Range& moab::LloydSmoother::elems ( ) const [inline]

Definition at line 87 of file LloydSmoother.hpp.

References myElems.

    {
        return myElems;
    }

Definition at line 94 of file LloydSmoother.hpp.

References fixedTag.

    {
        return fixedTag;
    }
void moab::LloydSmoother::fixed_tag ( Tag  fixed) [inline]

Definition at line 101 of file LloydSmoother.hpp.

References fixed, and fixedTag.

    {
        fixedTag = fixed;
    }

Definition at line 204 of file LloydSmoother.cpp.

References moab::ParallelComm::filter_pstatus(), moab::Skinner::find_skin(), fixed, fixedTag, iCreatedTag, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_OPAQUE, mbImpl, myElems, myPcomm, PSTATUS_GHOST, PSTATUS_INTERFACE, PSTATUS_NOT, moab::Range::size(), skin(), and smoab::UNION.

Referenced by perform_smooth().

{
    ErrorCode rval = MB_SUCCESS;

    // first, check for tag; if we don't have it, make one and mark skin as fixed
    if( !fixedTag )
    {
        unsigned char fixed = 0x0;
        rval = mbImpl->tag_get_handle( "", 1, MB_TYPE_OPAQUE, fixedTag, MB_TAG_DENSE | MB_TAG_CREAT, &fixed );MB_CHK_SET_ERR( rval, "Trouble making fixed tag" );
        iCreatedTag = true;

        // get the skin; get facets, because we might need to filter on shared entities
        Skinner skinner( mbImpl );
        Range skin, skin_verts;
        rval = skinner.find_skin( 0, myElems, false, skin );MB_CHK_SET_ERR( rval, "Unable to find skin" );

#ifdef MOAB_HAVE_MPI
        // need to do a little extra if we're working in parallel
        if( myPcomm )
        {
            // filter out ghost and interface facets
            rval = myPcomm->filter_pstatus( skin, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Failed to filter on shared status" );
        }
#endif
        // get the vertices from those entities
        rval = mbImpl->get_adjacencies( skin, 0, false, skin_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Trouble getting vertices" );

        // mark them fixed
        std::vector< unsigned char > marks( skin_verts.size(), 0x1 );
        rval = mbImpl->tag_set_data( fixedTag, skin_verts, &marks[0] );MB_CHK_SET_ERR( rval, "Unable to set tag on skin" );
    }

    return MB_SUCCESS;
}
Interface* moab::LloydSmoother::mb_impl ( ) [inline]

Definition at line 59 of file LloydSmoother.hpp.

References mbImpl.

    {
        return mbImpl;
    }
int moab::LloydSmoother::num_its ( ) [inline]

Definition at line 150 of file LloydSmoother.hpp.

References numIts.

Referenced by DeformMeshRemap::execute(), and main().

    {
        return numIts;
    }
void moab::LloydSmoother::num_its ( int  num) [inline]

Definition at line 154 of file LloydSmoother.hpp.

References numIts.

    {
        numIts = num;
    }

Definition at line 66 of file LloydSmoother.hpp.

References myPcomm.

    {
        return myPcomm;
    }
void moab::LloydSmoother::pcomm ( ParallelComm pc) [inline]

Definition at line 73 of file LloydSmoother.hpp.

References myPcomm.

    {
        myPcomm = pc;
    }

Definition at line 31 of file LloydSmoother.cpp.

References absTol, moab::Range::begin(), moab::ParallelComm::comm(), conn, coordsTag, moab::BoundBox::diagonal_length(), dim, moab::Range::empty(), moab::Range::end(), moab::ParallelComm::exchange_tags(), moab::ParallelComm::filter_pstatus(), fixedTag, initialize(), moab::Range::insert(), length(), moab::CN::MAX_NODES_PER_ELEMENT, MB_CHK_SET_ERR, MB_FAILURE, MB_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, mbImpl, myElems, myPcomm, numIts, PSTATUS_AND, PSTATUS_NOT, PSTATUS_NOT_OWNED, PSTATUS_SHARED, moab::ParallelComm::rank(), moab::Range::rbegin(), relTol, reportIts, moab::Range::size(), moab::ParallelComm::size(), smoab::UNION, and moab::BoundBox::update().

Referenced by DeformMeshRemap::execute(), main(), and run_local_smoother().

{
    ErrorCode rval;

    if( myElems.empty() ) { MB_SET_ERR( MB_FAILURE, "No elements specified to Lloyd smoother" ); }
    else if( mbImpl->dimension_from_handle( *myElems.begin() ) != mbImpl->dimension_from_handle( *myElems.rbegin() ) )
    {
        MB_SET_ERR( MB_FAILURE, "Elements of unequal dimension specified to Lloyd smoother" );
    }

    int dim = mbImpl->dimension_from_handle( *myElems.begin() );

    // first figure out tolerance to use
    if( 0 > absTol )
    {
        // no tolerance set - get one relative to bounding box around elements
        BoundBox bb;
        rval = bb.update( *mbImpl, myElems );MB_CHK_SET_ERR( rval, "Failed to compute bounding box around elements" );
        absTol = relTol * bb.diagonal_length();
    }

    // initialize if we need to
    rval = initialize();MB_CHK_SET_ERR( rval, "Failed to initialize" );

    // get all vertices
    Range verts;
    rval = mbImpl->get_adjacencies( myElems, 0, false, verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get all vertices" );

    // perform Lloyd relaxation:
    // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags

    // get all verts coords into tag; don't need to worry about filtering out fixed verts,
    // we'll just be setting to their fixed coords
    std::vector< double > vcentroids( 3 * verts.size() );
    if( !coordsTag )
    {
        rval = mbImpl->get_coords( verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get vert coords" );
    }
    else
    {
        rval = mbImpl->tag_get_data( coordsTag, verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get vert coords tag values" );
    }

    Tag centroid;
    rval = mbImpl->tag_get_handle( "", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_SET_ERR( rval, "Couldn't get tag handle" );
    rval = mbImpl->tag_set_data( centroid, verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed setting centroid tag" );

    Range owned_verts, shared_owned_verts;
#ifdef MOAB_HAVE_MPI
    // filter verts down to owned ones and get fixed tag for them
    if( myPcomm && myPcomm->size() > 1 )
    {
        rval = myPcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );MB_CHK_SET_ERR( rval, "Failed to filter on pstatus" );
        // get shared owned verts, for exchanging tags
        rval = myPcomm->filter_pstatus( owned_verts, PSTATUS_SHARED, PSTATUS_AND, -1, &shared_owned_verts );MB_CHK_SET_ERR( rval, "Failed to filter for shared owned" );
        // workaround: if no shared owned verts, put a non-shared one in the list, to prevent
        // exchanging tags for all shared entities
        if( shared_owned_verts.empty() ) shared_owned_verts.insert( *verts.begin() );
    }
    else
        owned_verts = verts;
#else
    owned_verts = verts;
#endif

    std::vector< unsigned char > fix_tag( owned_verts.size() );
    rval = mbImpl->tag_get_data( fixedTag, owned_verts, &fix_tag[0] );MB_CHK_SET_ERR( rval, "Failed to get fixed tag" );

    // now fill vcentroids array with positions of just owned vertices, since those are the ones
    // we're actually computing
    vcentroids.resize( 3 * owned_verts.size() );
    rval = mbImpl->tag_get_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get centroid tag" );

    // some declarations for later iterations
    std::vector< double > fcentroids( 3 * myElems.size() );  // fcentroids for element centroids
    std::vector< double > ctag(
        3 * CN::MAX_NODES_PER_ELEMENT );  // temporary coordinate storage for verts bounding an element
    const EntityHandle* conn;             // const ptr & size to elem connectivity
    int nconn;
    Range::iterator eit, vit;               // for iterating over elems, verts
    int e, v;                               // for indexing into centroid vectors
    std::vector< EntityHandle > adj_elems;  // used in vertex iteration

    // 2. while !converged
    double resid = DBL_MAX;
    numIts       = 0;
    while( resid > absTol )
    {
        numIts++;
        resid = 0.0;

        // 2a. foreach elem: centroid = sum(vertex centroids)/num_verts_in_cell
        for( eit = myElems.begin(), e = 0; eit != myElems.end(); ++eit, e++ )
        {
            // get verts for this elem
            rval = mbImpl->get_connectivity( *eit, conn, nconn );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
            // get centroid tags for those verts
            rval = mbImpl->tag_get_data( centroid, conn, nconn, &ctag[0] );MB_CHK_SET_ERR( rval, "Failed to get centroid" );
            fcentroids[3 * e + 0] = fcentroids[3 * e + 1] = fcentroids[3 * e + 2] = 0.0;
            for( v = 0; v < nconn; v++ )
            {
                fcentroids[3 * e + 0] += ctag[3 * v + 0];
                fcentroids[3 * e + 1] += ctag[3 * v + 1];
                fcentroids[3 * e + 2] += ctag[3 * v + 2];
            }
            for( v = 0; v < 3; v++ )
                fcentroids[3 * e + v] /= nconn;
        }
        rval = mbImpl->tag_set_data( centroid, myElems, &fcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set elem centroid" );

        // 2b. foreach owned vertex:
        for( vit = owned_verts.begin(), v = 0; vit != owned_verts.end(); ++vit, v++ )
        {
            // if !fixed
            if( fix_tag[v] ) continue;
            // vertex centroid = sum(cell centroids)/ncells
            adj_elems.clear();
            rval = mbImpl->get_adjacencies( &( *vit ), 1, dim, false, adj_elems );MB_CHK_SET_ERR( rval, "Failed getting adjs" );
            rval = mbImpl->tag_get_data( centroid, &adj_elems[0], adj_elems.size(), &fcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get elem centroid" );
            double vnew[] = { 0.0, 0.0, 0.0 };
            for( e = 0; e < (int)adj_elems.size(); e++ )
            {
                vnew[0] += fcentroids[3 * e + 0];
                vnew[1] += fcentroids[3 * e + 1];
                vnew[2] += fcentroids[3 * e + 2];
            }
            for( e = 0; e < 3; e++ )
                vnew[e] /= adj_elems.size();
            double delta = ( CartVect( vnew ) - CartVect( &vcentroids[3 * v] ) ).length();
            resid        = ( v ? std::max( resid, delta ) : delta );
            for( e = 0; e < 3; e++ )
                vcentroids[3 * v + e] = vnew[e];
        }

        // set the centroid tag; having them only in vcentroids array isn't enough, as vertex
        // centroids are accessed randomly in loop over faces
        rval = mbImpl->tag_set_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set vertex centroid" );

#ifdef MOAB_HAVE_MPI
        // 2c. exchange tags on owned verts
        if( myPcomm && myPcomm->size() > 1 )
        {
            rval = myPcomm->exchange_tags( centroid, shared_owned_verts );MB_CHK_SET_ERR( rval, "Failed to exchange tags" );
        }
#endif

        if( reportIts && !( numIts % reportIts ) )
        {
            double global_max = resid;
#ifdef MOAB_HAVE_MPI
            // global reduce for maximum delta, then report it
            if( myPcomm && myPcomm->size() > 1 )
                MPI_Reduce( &resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm() );
            if( !myPcomm || !myPcomm->rank() )
#endif
                std::cout << "Max residual = " << global_max << std::endl;
        }

    }  // end while

    // write the tag back onto vertex coordinates
    if( !coordsTag )
    {
        rval = mbImpl->set_coords( owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set vertex coords" );
    }
    else
    {
        rval = mbImpl->tag_set_data( coordsTag, owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set vert coords tag values" );
    }

    return MB_SUCCESS;
}
double moab::LloydSmoother::rel_tol ( ) [inline]

Definition at line 136 of file LloydSmoother.hpp.

References relTol.

    {
        return relTol;
    }
void moab::LloydSmoother::rel_tol ( double  tol) [inline]

Definition at line 143 of file LloydSmoother.hpp.

References relTol.

    {
        relTol = tol;
    }

Definition at line 161 of file LloydSmoother.hpp.

References reportIts.

Referenced by main().

    {
        return reportIts;
    }
void moab::LloydSmoother::report_its ( int  num) [inline]

Definition at line 165 of file LloydSmoother.hpp.

References reportIts.

    {
        reportIts = num;
    }

Member Data Documentation

double moab::LloydSmoother::absTol [private]

Definition at line 190 of file LloydSmoother.hpp.

Referenced by abs_tol(), and perform_smooth().

Definition at line 184 of file LloydSmoother.hpp.

Referenced by coords_tag(), and perform_smooth().

Definition at line 199 of file LloydSmoother.hpp.

Referenced by initialize(), and ~LloydSmoother().

Interface* moab::LloydSmoother::mbImpl [private]

Definition at line 175 of file LloydSmoother.hpp.

Referenced by initialize(), mb_impl(), perform_smooth(), and ~LloydSmoother().

Definition at line 181 of file LloydSmoother.hpp.

Referenced by elems(), initialize(), and perform_smooth().

Definition at line 178 of file LloydSmoother.hpp.

Referenced by initialize(), pcomm(), and perform_smooth().

Definition at line 196 of file LloydSmoother.hpp.

Referenced by num_its(), and perform_smooth().

double moab::LloydSmoother::relTol [private]

Definition at line 190 of file LloydSmoother.hpp.

Referenced by perform_smooth(), and rel_tol().

Definition at line 193 of file LloydSmoother.hpp.

Referenced by perform_smooth(), and report_its().

List of all members.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines