Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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 ()
Interfacemb_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

InterfacembImpl
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,
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 ErrorCode, fixedTag, iCreatedTag, MB_CHK_SET_ERR_RET, mbImpl, and moab::Interface::tag_delete().

{
    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 127 of file LloydSmoother.hpp.

References absTol.

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

Definition at line 134 of file LloydSmoother.hpp.

References absTol.

    {
        absTol = tol;
    }

Definition at line 113 of file LloydSmoother.hpp.

References coordsTag.

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

Definition at line 120 of file LloydSmoother.hpp.

References coordsTag.

    {
        coordsTag = coords;
    }

Definition at line 85 of file LloydSmoother.hpp.

References myElems.

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

Definition at line 92 of file LloydSmoother.hpp.

References myElems.

    {
        return myElems;
    }

Definition at line 99 of file LloydSmoother.hpp.

References fixedTag.

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

Definition at line 106 of file LloydSmoother.hpp.

References fixedTag.

    {
        fixedTag = fixed;
    }

Definition at line 207 of file LloydSmoother.cpp.

References ErrorCode, moab::ParallelComm::filter_pstatus(), moab::Skinner::find_skin(), fixedTag, moab::Interface::get_adjacencies(), 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(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), and moab::Interface::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;
}

Definition at line 64 of file LloydSmoother.hpp.

References mbImpl.

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

Definition at line 155 of file LloydSmoother.hpp.

References numIts.

Referenced by DeformMeshRemap::execute().

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

Definition at line 159 of file LloydSmoother.hpp.

References numIts.

    {
        numIts = num;
    }

Definition at line 71 of file LloydSmoother.hpp.

References myPcomm.

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

Definition at line 78 of file LloydSmoother.hpp.

References myPcomm.

    {
        myPcomm = pc;
    }

Definition at line 31 of file LloydSmoother.cpp.

References absTol, moab::Range::begin(), moab::ParallelComm::comm(), coordsTag, moab::BoundBox::diagonal_length(), dim, moab::Interface::dimension_from_handle(), moab::Range::empty(), moab::Range::end(), ErrorCode, moab::ParallelComm::exchange_tags(), moab::ParallelComm::filter_pstatus(), fixedTag, moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), initialize(), moab::Range::insert(), length(), moab::CN::MAX_NODES_PER_ELEMENT, MB_CHK_SET_ERR, 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::Interface::set_coords(), moab::Range::size(), moab::ParallelComm::size(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), moab::Interface::UNION, and moab::BoundBox::update().

Referenced by DeformMeshRemap::execute(), 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 141 of file LloydSmoother.hpp.

References relTol.

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

Definition at line 148 of file LloydSmoother.hpp.

References relTol.

    {
        relTol = tol;
    }

Definition at line 166 of file LloydSmoother.hpp.

References reportIts.

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

Definition at line 170 of file LloydSmoother.hpp.

References reportIts.

    {
        reportIts = num;
    }

Member Data Documentation

double moab::LloydSmoother::absTol [private]

Definition at line 195 of file LloydSmoother.hpp.

Referenced by abs_tol(), and perform_smooth().

Definition at line 189 of file LloydSmoother.hpp.

Referenced by coords_tag(), and perform_smooth().

Definition at line 204 of file LloydSmoother.hpp.

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

Definition at line 186 of file LloydSmoother.hpp.

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

Definition at line 183 of file LloydSmoother.hpp.

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

Definition at line 201 of file LloydSmoother.hpp.

Referenced by num_its(), and perform_smooth().

double moab::LloydSmoother::relTol [private]

Definition at line 195 of file LloydSmoother.hpp.

Referenced by perform_smooth(), and rel_tol().

Definition at line 198 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