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
|
