LCOV - code coverage report
Current view: top level - src - LloydSmoother.cpp (source / functions) Hit Total Coverage
Test: coverage_sk.info Lines: 93 100 93.0 %
Date: 2020-12-16 07:07:30 Functions: 6 6 100.0 %
Branches: 173 640 27.0 %

           Branch data     Line data    Source code
       1                 :            : #include "moab/LloydSmoother.hpp"
       2                 :            : #include "moab/Skinner.hpp"
       3                 :            : #include "moab/CN.hpp"
       4                 :            : #include "moab/CartVect.hpp"
       5                 :            : #include "moab/BoundBox.hpp"
       6                 :            : 
       7                 :            : #ifdef MOAB_HAVE_MPI
       8                 :            : #include "moab/ParallelComm.hpp"
       9                 :            : #include "MBParallelConventions.h"
      10                 :            : #endif
      11                 :            : 
      12                 :            : #include <iostream>
      13                 :            : 
      14                 :            : namespace moab
      15                 :            : {
      16                 :            : 
      17                 :          2 : LloydSmoother::LloydSmoother( Interface* impl, ParallelComm* pc, Range& elms, Tag ctag, Tag ftag, double at, double rt )
      18                 :            :     : mbImpl( impl ), myPcomm( pc ), myElems( elms ), coordsTag( ctag ), fixedTag( ftag ), absTol( at ), relTol( rt ),
      19                 :          2 :       reportIts( 0 ), numIts( 0 ), iCreatedTag( false )
      20                 :            : {
      21                 :          2 : }
      22                 :            : 
      23         [ +  - ]:          4 : LloydSmoother::~LloydSmoother()
      24                 :            : {
      25 [ +  - ][ +  - ]:          2 :     if( iCreatedTag && fixedTag )
      26                 :            :     {
      27 [ -  + ][ #  # ]:          2 :         ErrorCode rval = mbImpl->tag_delete( fixedTag );MB_CHK_SET_ERR_RET( rval, "Failed to delete the fixed tag" );
      28                 :            :     }
      29                 :          2 : }
      30                 :            : 
      31                 :          2 : ErrorCode LloydSmoother::perform_smooth()
      32                 :            : {
      33                 :            :     ErrorCode rval;
      34                 :            : 
      35 [ +  - ][ -  + ]:          2 :     if( myElems.empty() ) { MB_SET_ERR( MB_FAILURE, "No elements specified to Lloyd smoother" ); }
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      36 [ +  - ][ +  - ]:          2 :     else if( mbImpl->dimension_from_handle( *myElems.begin() ) != mbImpl->dimension_from_handle( *myElems.rbegin() ) )
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
                 [ -  + ]
      37                 :            :     {
      38 [ #  # ][ #  # ]:          0 :         MB_SET_ERR( MB_FAILURE, "Elements of unequal dimension specified to Lloyd smoother" );
         [ #  # ][ #  # ]
                 [ #  # ]
      39                 :            :     }
      40                 :            : 
      41 [ +  - ][ +  - ]:          2 :     int dim = mbImpl->dimension_from_handle( *myElems.begin() );
                 [ +  - ]
      42                 :            : 
      43                 :            :     // first figure out tolerance to use
      44         [ +  - ]:          2 :     if( 0 > absTol )
      45                 :            :     {
      46                 :            :         // no tolerance set - get one relative to bounding box around elements
      47         [ +  - ]:          2 :         BoundBox bb;
      48 [ +  - ][ -  + ]:          2 :         rval = bb.update( *mbImpl, myElems );MB_CHK_SET_ERR( rval, "Failed to compute bounding box around elements" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      49 [ +  - ][ +  - ]:          2 :         absTol = relTol * bb.diagonal_length();
      50                 :            :     }
      51                 :            : 
      52                 :            :     // initialize if we need to
      53 [ +  - ][ -  + ]:          2 :     rval = initialize();MB_CHK_SET_ERR( rval, "Failed to initialize" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      54                 :            : 
      55                 :            :     // get all vertices
      56         [ +  - ]:          2 :     Range verts;
      57 [ +  - ][ -  + ]:          2 :     rval = mbImpl->get_adjacencies( myElems, 0, false, verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get all vertices" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      58                 :            : 
      59                 :            :     // perform Lloyd relaxation:
      60                 :            :     // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
      61                 :            : 
      62                 :            :     // get all verts coords into tag; don't need to worry about filtering out fixed verts,
      63                 :            :     // we'll just be setting to their fixed coords
      64 [ +  - ][ +  - ]:          4 :     std::vector< double > vcentroids( 3 * verts.size() );
      65         [ +  + ]:          2 :     if( !coordsTag )
      66                 :            :     {
      67 [ +  - ][ +  - ]:          1 :         rval = mbImpl->get_coords( verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get vert coords" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      68                 :            :     }
      69                 :            :     else
      70                 :            :     {
      71 [ +  - ][ +  - ]:          1 :         rval = mbImpl->tag_get_data( coordsTag, verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get vert coords tag values" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      72                 :            :     }
      73                 :            : 
      74                 :            :     Tag centroid;
      75 [ +  - ][ -  + ]:          2 :     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" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      76 [ +  - ][ +  - ]:          2 :     rval = mbImpl->tag_set_data( centroid, verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed setting centroid tag" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      77                 :            : 
      78 [ +  - ][ +  - ]:          4 :     Range owned_verts, shared_owned_verts;
      79                 :            : #ifdef MOAB_HAVE_MPI
      80                 :            :     // filter verts down to owned ones and get fixed tag for them
      81 [ -  + ][ #  # ]:          2 :     if( myPcomm && myPcomm->size() > 1 )
         [ #  # ][ -  + ]
      82                 :            :     {
      83 [ #  # ][ #  # ]:          0 :         rval = myPcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );MB_CHK_SET_ERR( rval, "Failed to filter on pstatus" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      84                 :            :         // get shared owned verts, for exchanging tags
      85 [ #  # ][ #  # ]:          0 :         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" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      86                 :            :         // workaround: if no shared owned verts, put a non-shared one in the list, to prevent
      87                 :            :         // exchanging tags for all shared entities
      88 [ #  # ][ #  # ]:          0 :         if( shared_owned_verts.empty() ) shared_owned_verts.insert( *verts.begin() );
         [ #  # ][ #  # ]
                 [ #  # ]
      89                 :            :     }
      90                 :            :     else
      91         [ +  - ]:          2 :         owned_verts = verts;
      92                 :            : #else
      93                 :            :     owned_verts = verts;
      94                 :            : #endif
      95                 :            : 
      96 [ +  - ][ +  - ]:          4 :     std::vector< unsigned char > fix_tag( owned_verts.size() );
      97 [ +  - ][ +  - ]:          2 :     rval = mbImpl->tag_get_data( fixedTag, owned_verts, &fix_tag[0] );MB_CHK_SET_ERR( rval, "Failed to get fixed tag" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      98                 :            : 
      99                 :            :     // now fill vcentroids array with positions of just owned vertices, since those are the ones
     100                 :            :     // we're actually computing
     101 [ +  - ][ +  - ]:          2 :     vcentroids.resize( 3 * owned_verts.size() );
     102 [ +  - ][ +  - ]:          2 :     rval = mbImpl->tag_get_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get centroid tag" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     103                 :            : 
     104                 :            :     // some declarations for later iterations
     105 [ +  - ][ +  - ]:          4 :     std::vector< double > fcentroids( 3 * myElems.size() );  // fcentroids for element centroids
     106                 :            :     std::vector< double > ctag(
     107         [ +  - ]:          4 :         3 * CN::MAX_NODES_PER_ELEMENT );  // temporary coordinate storage for verts bounding an element
     108                 :            :     const EntityHandle* conn;             // const ptr & size to elem connectivity
     109                 :            :     int nconn;
     110 [ +  - ][ +  - ]:          2 :     Range::iterator eit, vit;               // for iterating over elems, verts
     111                 :            :     int e, v;                               // for indexing into centroid vectors
     112         [ +  - ]:          4 :     std::vector< EntityHandle > adj_elems;  // used in vertex iteration
     113                 :            : 
     114                 :            :     // 2. while !converged
     115                 :          2 :     double resid = DBL_MAX;
     116                 :          2 :     numIts       = 0;
     117         [ +  + ]:        285 :     while( resid > absTol )
     118                 :            :     {
     119                 :        283 :         numIts++;
     120                 :        283 :         resid = 0.0;
     121                 :            : 
     122                 :            :         // 2a. foreach elem: centroid = sum(vertex centroids)/num_verts_in_cell
     123 [ +  - ][ +  - ]:     174611 :         for( eit = myElems.begin(), e = 0; eit != myElems.end(); ++eit, e++ )
         [ +  - ][ +  - ]
                 [ +  + ]
     124                 :            :         {
     125                 :            :             // get verts for this elem
     126 [ +  - ][ +  - ]:     174328 :             rval = mbImpl->get_connectivity( *eit, conn, nconn );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     127                 :            :             // get centroid tags for those verts
     128 [ +  - ][ +  - ]:     174328 :             rval = mbImpl->tag_get_data( centroid, conn, nconn, &ctag[0] );MB_CHK_SET_ERR( rval, "Failed to get centroid" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     129 [ +  - ][ +  - ]:     174328 :             fcentroids[3 * e + 0] = fcentroids[3 * e + 1] = fcentroids[3 * e + 2] = 0.0;
                 [ +  - ]
     130         [ +  + ]:     697312 :             for( v = 0; v < nconn; v++ )
     131                 :            :             {
     132 [ +  - ][ +  - ]:     522984 :                 fcentroids[3 * e + 0] += ctag[3 * v + 0];
     133 [ +  - ][ +  - ]:     522984 :                 fcentroids[3 * e + 1] += ctag[3 * v + 1];
     134 [ +  - ][ +  - ]:     522984 :                 fcentroids[3 * e + 2] += ctag[3 * v + 2];
     135                 :            :             }
     136         [ +  + ]:     697312 :             for( v = 0; v < 3; v++ )
     137         [ +  - ]:     522984 :                 fcentroids[3 * e + v] /= nconn;
     138                 :            :         }
     139 [ +  - ][ +  - ]:        283 :         rval = mbImpl->tag_set_data( centroid, myElems, &fcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set elem centroid" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     140                 :            : 
     141                 :            :         // 2b. foreach owned vertex:
     142 [ +  - ][ +  - ]:      97635 :         for( vit = owned_verts.begin(), v = 0; vit != owned_verts.end(); ++vit, v++ )
         [ +  - ][ +  - ]
                 [ +  + ]
     143                 :            :         {
     144                 :            :             // if !fixed
     145 [ +  - ][ +  + ]:      97352 :             if( fix_tag[v] ) continue;
     146                 :            :             // vertex centroid = sum(cell centroids)/ncells
     147                 :      76976 :             adj_elems.clear();
     148 [ +  - ][ +  - ]:      76976 :             rval = mbImpl->get_adjacencies( &( *vit ), 1, dim, false, adj_elems );MB_CHK_SET_ERR( rval, "Failed getting adjs" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     149 [ +  - ][ +  - ]:      76976 :             rval = mbImpl->tag_get_data( centroid, &adj_elems[0], adj_elems.size(), &fcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to get elem centroid" );
         [ +  - ][ -  + ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     150                 :      76976 :             double vnew[] = { 0.0, 0.0, 0.0 };
     151         [ +  + ]:     541662 :             for( e = 0; e < (int)adj_elems.size(); e++ )
     152                 :            :             {
     153         [ +  - ]:     464686 :                 vnew[0] += fcentroids[3 * e + 0];
     154         [ +  - ]:     464686 :                 vnew[1] += fcentroids[3 * e + 1];
     155         [ +  - ]:     464686 :                 vnew[2] += fcentroids[3 * e + 2];
     156                 :            :             }
     157         [ +  + ]:     307904 :             for( e = 0; e < 3; e++ )
     158                 :     230928 :                 vnew[e] /= adj_elems.size();
     159 [ +  - ][ +  - ]:      76976 :             double delta = ( CartVect( vnew ) - CartVect( &vcentroids[3 * v] ) ).length();
         [ +  - ][ +  - ]
                 [ +  - ]
     160 [ +  - ][ +  - ]:      76976 :             resid        = ( v ? std::max( resid, delta ) : delta );
     161         [ +  + ]:     307904 :             for( e = 0; e < 3; e++ )
     162         [ +  - ]:     230928 :                 vcentroids[3 * v + e] = vnew[e];
     163                 :            :         }
     164                 :            : 
     165                 :            :         // set the centroid tag; having them only in vcentroids array isn't enough, as vertex
     166                 :            :         // centroids are accessed randomly in loop over faces
     167 [ +  - ][ +  - ]:        283 :         rval = mbImpl->tag_set_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set vertex centroid" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     168                 :            : 
     169                 :            : #ifdef MOAB_HAVE_MPI
     170                 :            :         // 2c. exchange tags on owned verts
     171 [ -  + ][ #  # ]:        283 :         if( myPcomm && myPcomm->size() > 1 )
         [ #  # ][ -  + ]
     172                 :            :         {
     173 [ #  # ][ #  # ]:          0 :             rval = myPcomm->exchange_tags( centroid, shared_owned_verts );MB_CHK_SET_ERR( rval, "Failed to exchange tags" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     174                 :            :         }
     175                 :            : #endif
     176                 :            : 
     177 [ +  - ][ +  + ]:        283 :         if( reportIts && !( numIts % reportIts ) )
     178                 :            :         {
     179                 :         28 :             double global_max = resid;
     180                 :            : #ifdef MOAB_HAVE_MPI
     181                 :            :             // global reduce for maximum delta, then report it
     182 [ -  + ][ #  # ]:         28 :             if( myPcomm && myPcomm->size() > 1 )
         [ #  # ][ -  + ]
     183 [ #  # ][ #  # ]:          0 :                 MPI_Reduce( &resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm() );
     184 [ -  + ][ #  # ]:         28 :             if( !myPcomm || !myPcomm->rank() )
         [ #  # ][ +  - ]
     185                 :            : #endif
     186 [ +  - ][ +  - ]:         28 :                 std::cout << "Max residual = " << global_max << std::endl;
                 [ +  - ]
     187                 :            :         }
     188                 :            : 
     189                 :            :     }  // end while
     190                 :            : 
     191                 :            :     // write the tag back onto vertex coordinates
     192         [ +  + ]:          2 :     if( !coordsTag )
     193                 :            :     {
     194 [ +  - ][ +  - ]:          1 :         rval = mbImpl->set_coords( owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set vertex coords" );
         [ +  - ][ -  + ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     195                 :            :     }
     196                 :            :     else
     197                 :            :     {
     198 [ +  - ][ +  - ]:          1 :         rval = mbImpl->tag_set_data( coordsTag, owned_verts, &vcentroids[0] );MB_CHK_SET_ERR( rval, "Failed to set vert coords tag values" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     199                 :            :     }
     200                 :            : 
     201                 :          4 :     return MB_SUCCESS;
     202                 :            : }
     203                 :            : 
     204                 :          2 : ErrorCode LloydSmoother::initialize()
     205                 :            : {
     206                 :          2 :     ErrorCode rval = MB_SUCCESS;
     207                 :            : 
     208                 :            :     // first, check for tag; if we don't have it, make one and mark skin as fixed
     209         [ +  - ]:          2 :     if( !fixedTag )
     210                 :            :     {
     211                 :          2 :         unsigned char fixed = 0x0;
     212 [ +  - ][ -  + ]:          2 :         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" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     213                 :          2 :         iCreatedTag = true;
     214                 :            : 
     215                 :            :         // get the skin; get facets, because we might need to filter on shared entities
     216         [ +  - ]:          2 :         Skinner skinner( mbImpl );
     217 [ +  - ][ +  - ]:          4 :         Range skin, skin_verts;
         [ +  - ][ +  - ]
     218 [ +  - ][ -  + ]:          2 :         rval = skinner.find_skin( 0, myElems, false, skin );MB_CHK_SET_ERR( rval, "Unable to find skin" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     219                 :            : 
     220                 :            : #ifdef MOAB_HAVE_MPI
     221                 :            :         // need to do a little extra if we're working in parallel
     222         [ -  + ]:          2 :         if( myPcomm )
     223                 :            :         {
     224                 :            :             // filter out ghost and interface facets
     225 [ #  # ][ #  # ]:          0 :             rval = myPcomm->filter_pstatus( skin, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Failed to filter on shared status" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     226                 :            :         }
     227                 :            : #endif
     228                 :            :         // get the vertices from those entities
     229 [ +  - ][ -  + ]:          2 :         rval = mbImpl->get_adjacencies( skin, 0, false, skin_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Trouble getting vertices" );
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     230                 :            : 
     231                 :            :         // mark them fixed
     232 [ +  - ][ +  - ]:          4 :         std::vector< unsigned char > marks( skin_verts.size(), 0x1 );
                 [ +  - ]
     233 [ +  - ][ +  - ]:          4 :         rval = mbImpl->tag_set_data( fixedTag, skin_verts, &marks[0] );MB_CHK_SET_ERR( rval, "Unable to set tag on skin" );
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ +  - ]
     234                 :            :     }
     235                 :            : 
     236                 :          2 :     return MB_SUCCESS;
     237                 :            : }
     238                 :            : 
     239 [ +  - ][ +  - ]:          4 : }  // namespace moab

Generated by: LCOV version 1.11