#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

Detailed Description

Examples:: DeformMeshRemap.cpp.

Definition at line 25 of file LloydSmoother.hpp.

Constructor & Destructor Documentation

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

moab::LloydSmoother::~LloydSmoother ( )

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

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

ErrorCode moab::LloydSmoother::perform_smooth ( )

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