1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
#include "moab/LloydSmoother.hpp"
#include "moab/Skinner.hpp"
#include "moab/CN.hpp"
#include "moab/CartVect.hpp"
#include "moab/BoundBox.hpp"

#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#endif

#include <iostream>

namespace moab
{

LloydSmoother::LloydSmoother( Interface* impl, ParallelComm* pc, Range& elms, Tag ctag, Tag ftag, double at, double rt )
    : mbImpl( impl ), myPcomm( pc ), myElems( elms ), coordsTag( ctag ), fixedTag( ftag ), absTol( at ), relTol( rt ),
      reportIts( 0 ), numIts( 0 ), iCreatedTag( false )
{
}

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

ErrorCode LloydSmoother::perform_smooth()<--- The function 'perform_smooth' is never used.
{
    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;
}

ErrorCode LloydSmoother::initialize()
{
    ErrorCode rval = MB_SUCCESS;<--- Variable 'rval' is assigned a value that is never used.

    // 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;
}

}  // namespace moab