MOAB: Mesh Oriented datABase  (version 5.2.1)
LloydSmoother.cpp
Go to the documentation of this file.
00001 #include "moab/LloydSmoother.hpp"
00002 #include "moab/Skinner.hpp"
00003 #include "moab/CN.hpp"
00004 #include "moab/CartVect.hpp"
00005 #include "moab/BoundBox.hpp"
00006 
00007 #ifdef MOAB_HAVE_MPI
00008 #include "moab/ParallelComm.hpp"
00009 #include "MBParallelConventions.h"
00010 #endif
00011 
00012 #include <iostream>
00013 
00014 namespace moab
00015 {
00016 
00017 LloydSmoother::LloydSmoother( Interface* impl, ParallelComm* pc, Range& elms, Tag ctag, Tag ftag, double at, double rt )
00018     : mbImpl( impl ), myPcomm( pc ), myElems( elms ), coordsTag( ctag ), fixedTag( ftag ), absTol( at ), relTol( rt ),
00019       reportIts( 0 ), numIts( 0 ), iCreatedTag( false )
00020 {
00021 }
00022 
00023 LloydSmoother::~LloydSmoother()
00024 {
00025     if( iCreatedTag && fixedTag )
00026     {
00027         ErrorCode rval = mbImpl->tag_delete( fixedTag );MB_CHK_SET_ERR_RET( rval, "Failed to delete the fixed tag" );
00028     }
00029 }
00030 
00031 ErrorCode LloydSmoother::perform_smooth()
00032 {
00033     ErrorCode rval;
00034 
00035     if( myElems.empty() ) { MB_SET_ERR( MB_FAILURE, "No elements specified to Lloyd smoother" ); }
00036     else if( mbImpl->dimension_from_handle( *myElems.begin() ) != mbImpl->dimension_from_handle( *myElems.rbegin() ) )
00037     {
00038         MB_SET_ERR( MB_FAILURE, "Elements of unequal dimension specified to Lloyd smoother" );
00039     }
00040 
00041     int dim = mbImpl->dimension_from_handle( *myElems.begin() );
00042 
00043     // first figure out tolerance to use
00044     if( 0 > absTol )
00045     {
00046         // no tolerance set - get one relative to bounding box around elements
00047         BoundBox bb;
00048         rval = bb.update( *mbImpl, myElems );MB_CHK_SET_ERR( rval, "Failed to compute bounding box around elements" );
00049         absTol = relTol * bb.diagonal_length();
00050     }
00051 
00052     // initialize if we need to
00053     rval = initialize();MB_CHK_SET_ERR( rval, "Failed to initialize" );
00054 
00055     // get all vertices
00056     Range verts;
00057     rval = mbImpl->get_adjacencies( myElems, 0, false, verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get all vertices" );
00058 
00059     // perform Lloyd relaxation:
00060     // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
00061 
00062     // get all verts coords into tag; don't need to worry about filtering out fixed verts,
00063     // we'll just be setting to their fixed coords
00064     std::vector< double > vcentroids( 3 * verts.size() );
00065     if( !coordsTag )
00066     {
00067         rval = mbImpl->get_coords( verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get vert coords" );
00068     }
00069     else
00070     {
00071         rval = mbImpl->tag_get_data( coordsTag, verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get vert coords tag values" );
00072     }
00073 
00074     Tag centroid;
00075     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" );
00076     rval = mbImpl->tag_set_data( centroid, verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed setting centroid tag" );
00077 
00078     Range owned_verts, shared_owned_verts;
00079 #ifdef MOAB_HAVE_MPI
00080     // filter verts down to owned ones and get fixed tag for them
00081     if( myPcomm && myPcomm->size() > 1 )
00082     {
00083         rval = myPcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );MB_CHK_SET_ERR( rval, "Failed to filter on pstatus" );
00084         // get shared owned verts, for exchanging tags
00085         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" );
00086         // workaround: if no shared owned verts, put a non-shared one in the list, to prevent
00087         // exchanging tags for all shared entities
00088         if( shared_owned_verts.empty() ) shared_owned_verts.insert( *verts.begin() );
00089     }
00090     else
00091         owned_verts = verts;
00092 #else
00093     owned_verts = verts;
00094 #endif
00095 
00096     std::vector< unsigned char > fix_tag( owned_verts.size() );
00097     rval = mbImpl->tag_get_data( fixedTag, owned_verts, &fix_tag[0] );MB_CHK_SET_ERR( rval, "Failed to get fixed tag" );
00098 
00099     // now fill vcentroids array with positions of just owned vertices, since those are the ones
00100     // we're actually computing
00101     vcentroids.resize( 3 * owned_verts.size() );
00102     rval = mbImpl->tag_get_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get centroid tag" );
00103 
00104     // some declarations for later iterations
00105     std::vector< double > fcentroids( 3 * myElems.size() );  // fcentroids for element centroids
00106     std::vector< double > ctag(
00107         3 * CN::MAX_NODES_PER_ELEMENT );  // temporary coordinate storage for verts bounding an element
00108     const EntityHandle* conn;             // const ptr & size to elem connectivity
00109     int nconn;
00110     Range::iterator eit, vit;               // for iterating over elems, verts
00111     int e, v;                               // for indexing into centroid vectors
00112     std::vector< EntityHandle > adj_elems;  // used in vertex iteration
00113 
00114     // 2. while !converged
00115     double resid = DBL_MAX;
00116     numIts       = 0;
00117     while( resid > absTol )
00118     {
00119         numIts++;
00120         resid = 0.0;
00121 
00122         // 2a. foreach elem: centroid = sum(vertex centroids)/num_verts_in_cell
00123         for( eit = myElems.begin(), e = 0; eit != myElems.end(); ++eit, e++ )
00124         {
00125             // get verts for this elem
00126             rval = mbImpl->get_connectivity( *eit, conn, nconn );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
00127             // get centroid tags for those verts
00128             rval = mbImpl->tag_get_data( centroid, conn, nconn, &ctag[0] );MB_CHK_SET_ERR( rval, "Failed to get centroid" );
00129             fcentroids[3 * e + 0] = fcentroids[3 * e + 1] = fcentroids[3 * e + 2] = 0.0;
00130             for( v = 0; v < nconn; v++ )
00131             {
00132                 fcentroids[3 * e + 0] += ctag[3 * v + 0];
00133                 fcentroids[3 * e + 1] += ctag[3 * v + 1];
00134                 fcentroids[3 * e + 2] += ctag[3 * v + 2];
00135             }
00136             for( v = 0; v < 3; v++ )
00137                 fcentroids[3 * e + v] /= nconn;
00138         }
00139         rval = mbImpl->tag_set_data( centroid, myElems, &fcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set elem centroid" );
00140 
00141         // 2b. foreach owned vertex:
00142         for( vit = owned_verts.begin(), v = 0; vit != owned_verts.end(); ++vit, v++ )
00143         {
00144             // if !fixed
00145             if( fix_tag[v] ) continue;
00146             // vertex centroid = sum(cell centroids)/ncells
00147             adj_elems.clear();
00148             rval = mbImpl->get_adjacencies( &( *vit ), 1, dim, false, adj_elems );MB_CHK_SET_ERR( rval, "Failed getting adjs" );
00149             rval = mbImpl->tag_get_data( centroid, &adj_elems[0], adj_elems.size(), &fcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get elem centroid" );
00150             double vnew[] = { 0.0, 0.0, 0.0 };
00151             for( e = 0; e < (int)adj_elems.size(); e++ )
00152             {
00153                 vnew[0] += fcentroids[3 * e + 0];
00154                 vnew[1] += fcentroids[3 * e + 1];
00155                 vnew[2] += fcentroids[3 * e + 2];
00156             }
00157             for( e = 0; e < 3; e++ )
00158                 vnew[e] /= adj_elems.size();
00159             double delta = ( CartVect( vnew ) - CartVect( &vcentroids[3 * v] ) ).length();
00160             resid        = ( v ? std::max( resid, delta ) : delta );
00161             for( e = 0; e < 3; e++ )
00162                 vcentroids[3 * v + e] = vnew[e];
00163         }
00164 
00165         // set the centroid tag; having them only in vcentroids array isn't enough, as vertex
00166         // centroids are accessed randomly in loop over faces
00167         rval = mbImpl->tag_set_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set vertex centroid" );
00168 
00169 #ifdef MOAB_HAVE_MPI
00170         // 2c. exchange tags on owned verts
00171         if( myPcomm && myPcomm->size() > 1 )
00172         {
00173             rval = myPcomm->exchange_tags( centroid, shared_owned_verts );MB_CHK_SET_ERR( rval, "Failed to exchange tags" );
00174         }
00175 #endif
00176 
00177         if( reportIts && !( numIts % reportIts ) )
00178         {
00179             double global_max = resid;
00180 #ifdef MOAB_HAVE_MPI
00181             // global reduce for maximum delta, then report it
00182             if( myPcomm && myPcomm->size() > 1 )
00183                 MPI_Reduce( &resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm() );
00184             if( !myPcomm || !myPcomm->rank() )
00185 #endif
00186                 std::cout << "Max residual = " << global_max << std::endl;
00187         }
00188 
00189     }  // end while
00190 
00191     // write the tag back onto vertex coordinates
00192     if( !coordsTag )
00193     {
00194         rval = mbImpl->set_coords( owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set vertex coords" );
00195     }
00196     else
00197     {
00198         rval = mbImpl->tag_set_data( coordsTag, owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set vert coords tag values" );
00199     }
00200 
00201     return MB_SUCCESS;
00202 }
00203 
00204 ErrorCode LloydSmoother::initialize()
00205 {
00206     ErrorCode rval = MB_SUCCESS;
00207 
00208     // first, check for tag; if we don't have it, make one and mark skin as fixed
00209     if( !fixedTag )
00210     {
00211         unsigned char fixed = 0x0;
00212         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" );
00213         iCreatedTag = true;
00214 
00215         // get the skin; get facets, because we might need to filter on shared entities
00216         Skinner skinner( mbImpl );
00217         Range skin, skin_verts;
00218         rval = skinner.find_skin( 0, myElems, false, skin );MB_CHK_SET_ERR( rval, "Unable to find skin" );
00219 
00220 #ifdef MOAB_HAVE_MPI
00221         // need to do a little extra if we're working in parallel
00222         if( myPcomm )
00223         {
00224             // filter out ghost and interface facets
00225             rval = myPcomm->filter_pstatus( skin, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Failed to filter on shared status" );
00226         }
00227 #endif
00228         // get the vertices from those entities
00229         rval = mbImpl->get_adjacencies( skin, 0, false, skin_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Trouble getting vertices" );
00230 
00231         // mark them fixed
00232         std::vector< unsigned char > marks( skin_verts.size(), 0x1 );
00233         rval = mbImpl->tag_set_data( fixedTag, skin_verts, &marks[0] );MB_CHK_SET_ERR( rval, "Unable to set tag on skin" );
00234     }
00235 
00236     return MB_SUCCESS;
00237 }
00238 
00239 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines