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