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