moab
|
#include <LloydSmoother.hpp>
Public Member Functions | |
LloydSmoother (Interface *impl) | |
LloydSmoother (Interface *impl, ParallelComm *pc, Range &elems, Tag cds_tag=0, Tag fixed_tag=0, double abs_tol=-1.0, double rel_tol=1.0e-6) | |
~LloydSmoother () | |
ErrorCode | perform_smooth () |
Interface * | mb_impl () |
ParallelComm * | pcomm () |
void | pcomm (ParallelComm *pc) |
Range & | elems () |
const Range & | elems () const |
Tag | fixed_tag () |
void | fixed_tag (Tag fixed) |
Tag | coords_tag () |
void | coords_tag (Tag coords) |
double | abs_tol () |
void | abs_tol (double tol) |
double | rel_tol () |
void | rel_tol (double tol) |
int | num_its () |
void | num_its (int num) |
int | report_its () |
void | report_its (int num) |
Private Member Functions | |
ErrorCode | initialize () |
Private Attributes | |
Interface * | mbImpl |
ParallelComm * | myPcomm |
Range | myElems |
Tag | coordsTag |
Tag | fixedTag |
double | absTol |
double | relTol |
Error * | errorHandler |
int | reportIts |
int | numIts |
bool | iCreatedTag |
Definition at line 25 of file LloydSmoother.hpp.
moab::LloydSmoother::LloydSmoother | ( | Interface * | impl | ) |
moab::LloydSmoother::LloydSmoother | ( | Interface * | impl, |
ParallelComm * | pc, | ||
Range & | elems, | ||
Tag | cds_tag = 0 , |
||
Tag | fixed_tag = 0 , |
||
double | abs_tol = -1.0 , |
||
double | rel_tol = 1.0e-6 |
||
) |
Definition at line 21 of file LloydSmoother.cpp.
: mbImpl(impl), myPcomm(pc), myElems(elms), coordsTag(ctag), fixedTag(ftag), absTol(at), relTol(rt), reportIts(0), numIts(0), iCreatedTag(false) { ErrorCode rval = mbImpl->query_interface(errorHandler); if (rval) throw rval; }
Definition at line 30 of file LloydSmoother.cpp.
{ if (iCreatedTag && fixedTag) { ErrorCode rval = mbImpl->tag_delete(fixedTag); if (rval) throw rval; } }
double moab::LloydSmoother::abs_tol | ( | ) | [inline] |
Definition at line 95 of file LloydSmoother.hpp.
{return absTol;}
void moab::LloydSmoother::abs_tol | ( | double | tol | ) | [inline] |
Definition at line 99 of file LloydSmoother.hpp.
{absTol = tol;}
Tag moab::LloydSmoother::coords_tag | ( | ) | [inline] |
Definition at line 87 of file LloydSmoother.hpp.
{return coordsTag;}
void moab::LloydSmoother::coords_tag | ( | Tag | coords | ) | [inline] |
Definition at line 91 of file LloydSmoother.hpp.
{coordsTag = coords;}
Range& moab::LloydSmoother::elems | ( | ) | [inline] |
Definition at line 71 of file LloydSmoother.hpp.
{return myElems;}
const Range& moab::LloydSmoother::elems | ( | ) | const [inline] |
Definition at line 75 of file LloydSmoother.hpp.
{return myElems;}
Tag moab::LloydSmoother::fixed_tag | ( | ) | [inline] |
Definition at line 79 of file LloydSmoother.hpp.
{return fixedTag;}
void moab::LloydSmoother::fixed_tag | ( | Tag | fixed | ) | [inline] |
Definition at line 83 of file LloydSmoother.hpp.
{fixedTag = fixed;}
ErrorCode moab::LloydSmoother::initialize | ( | ) | [private] |
Definition at line 205 of file LloydSmoother.cpp.
{ ErrorCode rval = MB_SUCCESS; // 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); RR("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); RR("Unable to find skin."); #ifdef USE_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); RR("Failed to filter on shared status."); } #endif // get the vertices from those entities rval = mbImpl->get_adjacencies(skin, 0, false, skin_verts, Interface::UNION); RR("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]); RR("Unable to set tag on skin."); } return MB_SUCCESS; }
Interface* moab::LloydSmoother::mb_impl | ( | ) | [inline] |
Definition at line 59 of file LloydSmoother.hpp.
{return mbImpl;}
int moab::LloydSmoother::num_its | ( | ) | [inline] |
Definition at line 111 of file LloydSmoother.hpp.
{return numIts;}
void moab::LloydSmoother::num_its | ( | int | num | ) | [inline] |
Definition at line 112 of file LloydSmoother.hpp.
{numIts = num;}
ParallelComm* moab::LloydSmoother::pcomm | ( | ) | [inline] |
Definition at line 63 of file LloydSmoother.hpp.
{return myPcomm;}
void moab::LloydSmoother::pcomm | ( | ParallelComm * | pc | ) | [inline] |
Definition at line 67 of file LloydSmoother.hpp.
{myPcomm = pc;}
Definition at line 38 of file LloydSmoother.cpp.
{ ErrorCode rval; if (myElems.empty()) { rval = MB_FAILURE; RR("No elements specified to Lloyd smoother."); } else if (mbImpl->dimension_from_handle(*myElems.begin()) != mbImpl->dimension_from_handle(*myElems.rbegin())) { rval = MB_FAILURE; RR("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); RR("Failed to compute bounding box around elements."); absTol = relTol * bb.diagonal_length(); } // initialize if we need to rval = initialize(); RR("Failed to initialize."); // get all vertices Range verts; rval = mbImpl->get_adjacencies(myElems, 0, false, verts, Interface::UNION); RR("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]); RR("Failed to get vert coords."); } else { rval = mbImpl->tag_get_data(coordsTag, verts, &vcentroids[0]); RR("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); RR("Couldn't get tag handle."); rval = mbImpl->tag_set_data(centroid, verts, &vcentroids[0]); RR("Failed setting centroid tag."); Range owned_verts, shared_owned_verts; #ifdef USE_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); RR("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); RR("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]); RR("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]); RR("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); RR("Failed to get connectivity."); // get centroid tags for those verts rval = mbImpl->tag_get_data(centroid, conn, nconn, &ctag[0]); RR("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]); RR("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); RR("Failed getting adjs."); rval = mbImpl->tag_get_data(centroid, &adj_elems[0], adj_elems.size(), &fcentroids[0]); RR("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]); RR("Failed to set vertex centroid."); #ifdef USE_MPI // 2c. exchange tags on owned verts if (myPcomm && myPcomm->size() > 1) { rval = myPcomm->exchange_tags(centroid, shared_owned_verts); RR("Failed to exchange tags."); } #endif if (reportIts && !(numIts%reportIts)) { double global_max = resid; #ifdef USE_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]); RR("Failed to set vertex coords."); } else { rval = mbImpl->tag_set_data(coordsTag, owned_verts, &vcentroids[0]); RR("Failed to set vert coords tag values."); } return MB_SUCCESS; }
double moab::LloydSmoother::rel_tol | ( | ) | [inline] |
Definition at line 103 of file LloydSmoother.hpp.
{return relTol;}
void moab::LloydSmoother::rel_tol | ( | double | tol | ) | [inline] |
Definition at line 107 of file LloydSmoother.hpp.
{relTol = tol;}
int moab::LloydSmoother::report_its | ( | ) | [inline] |
Definition at line 116 of file LloydSmoother.hpp.
{return reportIts;}
void moab::LloydSmoother::report_its | ( | int | num | ) | [inline] |
Definition at line 117 of file LloydSmoother.hpp.
{reportIts = num;}
double moab::LloydSmoother::absTol [private] |
Definition at line 141 of file LloydSmoother.hpp.
Tag moab::LloydSmoother::coordsTag [private] |
Definition at line 135 of file LloydSmoother.hpp.
Error* moab::LloydSmoother::errorHandler [private] |
Definition at line 144 of file LloydSmoother.hpp.
Tag moab::LloydSmoother::fixedTag [private] |
Definition at line 138 of file LloydSmoother.hpp.
bool moab::LloydSmoother::iCreatedTag [private] |
Definition at line 153 of file LloydSmoother.hpp.
Interface* moab::LloydSmoother::mbImpl [private] |
Definition at line 126 of file LloydSmoother.hpp.
Range moab::LloydSmoother::myElems [private] |
Definition at line 132 of file LloydSmoother.hpp.
ParallelComm* moab::LloydSmoother::myPcomm [private] |
Definition at line 129 of file LloydSmoother.hpp.
int moab::LloydSmoother::numIts [private] |
Definition at line 150 of file LloydSmoother.hpp.
double moab::LloydSmoother::relTol [private] |
Definition at line 141 of file LloydSmoother.hpp.
int moab::LloydSmoother::reportIts [private] |
Definition at line 147 of file LloydSmoother.hpp.