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