moab
moab::LloydSmoother Class Reference

#include <LloydSmoother.hpp>

List of all members.

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 ()
Interfacemb_impl ()
ParallelCommpcomm ()
void pcomm (ParallelComm *pc)
Rangeelems ()
const Rangeelems () 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

InterfacembImpl
ParallelCommmyPcomm
Range myElems
Tag coordsTag
Tag fixedTag
double absTol
double relTol
ErrorerrorHandler
int reportIts
int numIts
bool iCreatedTag

Detailed Description

Examples:
DeformMeshRemap.cpp.

Definition at line 25 of file LloydSmoother.hpp.


Constructor & Destructor Documentation

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;
  }
}

Member Function Documentation

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;}

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;}

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;}

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;}

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;
}

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;}

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;}

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;}

Member Data Documentation

double moab::LloydSmoother::absTol [private]

Definition at line 141 of file LloydSmoother.hpp.

Definition at line 135 of file LloydSmoother.hpp.

Definition at line 144 of file LloydSmoother.hpp.

Definition at line 138 of file LloydSmoother.hpp.

Definition at line 153 of file LloydSmoother.hpp.

Definition at line 126 of file LloydSmoother.hpp.

Definition at line 132 of file LloydSmoother.hpp.

Definition at line 129 of file LloydSmoother.hpp.

Definition at line 150 of file LloydSmoother.hpp.

double moab::LloydSmoother::relTol [private]

Definition at line 141 of file LloydSmoother.hpp.

Definition at line 147 of file LloydSmoother.hpp.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines