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