00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceAToCPointLocator.hpp"
00032 #include "SundanceOut.hpp"
00033 #include "PlayaTabs.hpp"
00034 #include "SundanceGeomUtils.hpp"
00035 #include "Teuchos_ScalarTraits.hpp"
00036 #include <queue>
00037
00038 using namespace Sundance;
00039 using namespace Sundance;
00040 using namespace Sundance;
00041 using namespace Sundance;
00042 using namespace Sundance;
00043 using namespace Sundance;
00044 using namespace Sundance;
00045 using namespace Teuchos;
00046
00047
00048
00049
00050
00051 static Time& pointLocatorCtorTimer()
00052 {
00053 static RCP<Time> rtn
00054 = TimeMonitor::getNewTimer("point locator ctor");
00055 return *rtn;
00056 }
00057
00058 AToCPointLocator::AToCPointLocator(const Mesh& mesh,
00059 const CellFilter& subdomain,
00060 const std::vector<int>& nx)
00061 : dim_(mesh.spatialDim()),
00062 mesh_(mesh),
00063 nFacets_(mesh.numFacets(dim_, 0, 0)),
00064 nx_(nx),
00065 low_(nx.size(), 1.0/ScalarTraits<double>::sfmin()),
00066 high_(nx.size(), -1.0/ScalarTraits<double>::sfmin()),
00067 dx_(nx.size()),
00068 table_(),
00069 subdomain_(subdomain),
00070 neighborSet_()
00071 {
00072 TimeMonitor timer(pointLocatorCtorTimer());
00073
00074
00075 neighborSet_.resize(mesh.numCells(dim_));
00076
00077
00078 CellSet cells = subdomain.getCells(mesh);
00079
00080 for (CellIterator i = cells.begin(); i!= cells.end(); i++)
00081 {
00082 int cellLID = *i;
00083 Array<int> facetLIDs;
00084 Array<int> facetOri;
00085 mesh.getFacetArray(dim_, cellLID, 0, facetLIDs, facetOri);
00086 for (int f=0; f<facetLIDs.size(); f++)
00087 {
00088 Point x = mesh.nodePosition(facetLIDs[f]);
00089 for (int d=0; d<dim_; d++)
00090 {
00091 if (x[d] < low_[d]) low_[d] = x[d];
00092 if (x[d] > high_[d]) high_[d] = x[d];
00093 }
00094 }
00095 }
00096
00097
00098 for (int d=0; d<dim_; d++)
00099 {
00100 low_[d] -= 0.01 * (high_[d] - low_[d]);
00101 high_[d] += 0.01 * (high_[d] - low_[d]);
00102 }
00103
00104
00105
00106 int s = 1;
00107 for (int d=0; d<dim_; d++)
00108 {
00109 dx_[d] = (high_[d] - low_[d])/nx_[d];
00110 s *= nx_[d];
00111 }
00112
00113
00114 table_ = rcp(new Array<int>(s, -1));
00115
00116
00117 Array<int> lowIndex;
00118 Array<int> highIndex;
00119 for (CellIterator i = cells.begin(); i!= cells.end(); i++)
00120 {
00121 int cellLID = *i;
00122 getGridRange(mesh, dim_, cellLID, lowIndex, highIndex);
00123 if (dim_==2)
00124 {
00125 for (int ix=lowIndex[0]; ix<=highIndex[0]; ix++)
00126 {
00127 for (int iy=lowIndex[1]; iy<=highIndex[1]; iy++)
00128 {
00129 int index = nx_[1]*ix + iy;
00130 (*table_)[index] = cellLID;
00131 }
00132 }
00133 }
00134 else
00135 {
00136 TEUCHOS_TEST_FOR_EXCEPT(true);
00137 }
00138 }
00139 }
00140
00141 int AToCPointLocator::getGridIndex(const double* x) const
00142 {
00143 int index = 0;
00144 for (int d=0; d<dim_; d++)
00145 {
00146 double r = (x[d] - low_[d])/dx_[d];
00147 int ix = (int) floor(r);
00148 index = index*nx_[d] + ix;
00149 }
00150
00151 return index;
00152 }
00153
00154
00155 void AToCPointLocator::getGridRange(const Mesh& mesh, int cellDim, int cellLID,
00156 Array<int>& lowIndex, Array<int>& highIndex) const
00157 {
00158 Array<int> facetLIDs;
00159 Array<int> facetOri;
00160 mesh.getFacetArray(cellDim, cellLID, 0, facetLIDs, facetOri);
00161
00162 lowIndex.resize(cellDim);
00163 highIndex.resize(cellDim);
00164 for (int d=0; d<cellDim; d++)
00165 {
00166 highIndex[d] = -1;
00167 lowIndex[d] = 1000000;
00168
00169 }
00170
00171 for (int f=0; f<facetLIDs.size(); f++)
00172 {
00173 Point x = mesh.nodePosition(facetLIDs[f]);
00174 for (int d=0; d<cellDim; d++)
00175 {
00176 double r = (x[d] - low_[d])/dx_[d];
00177 int ix = (int) floor(r);
00178 if (ix < lowIndex[d]) lowIndex[d] = ix;
00179 if (ix > highIndex[d]) highIndex[d] = ix;
00180 }
00181 }
00182 }
00183
00184
00185 void AToCPointLocator::fillMaximalNeighbors(int cellLID,
00186 const int* facetLID) const
00187 {
00188 if (neighborSet_[cellLID].get()==0)
00189 {
00190 neighborSet_[cellLID] = rcp(new Set<int>());
00191
00192 for (int f=0; f<nFacets_; f++)
00193 {
00194 Array<int> cofacets;
00195 mesh_.getCofacets(0, facetLID[f], dim_, cofacets);
00196 for (int c=0; c<cofacets.size(); c++)
00197 {
00198 if (cofacets[c] != cellLID) neighborSet_[cellLID]->put(cofacets[c]);
00199 }
00200 }
00201 }
00202 }
00203
00204
00205 bool AToCPointLocator::cellContainsPoint(int cellLID,
00206 const double* x,
00207 const int* facetLID) const
00208 {
00209 if (dim_==2)
00210 {
00211 const double* A = mesh_.nodePositionView(facetLID[0]);
00212 const double* B = mesh_.nodePositionView(facetLID[1]);
00213 const double* C = mesh_.nodePositionView(facetLID[2]);
00214
00215
00216 double sign = orient2D(A, B, C);
00217 if (sign > 0.0)
00218 {
00219 double s1 = orient2D(A, B, x);
00220 double s2 = orient2D(B, C, x);
00221 double s3 = orient2D(C, A, x);
00222 if (s1 >= 0.0 && s2 >= 0.0 && s3 >= 0.0) return true;
00223 return false;
00224 }
00225 else
00226 {
00227 double s1 = orient2D(A, C, x);
00228 double s2 = orient2D(C, B, x);
00229 double s3 = orient2D(B, A, x);
00230 if (s1 >= 0.0 && s2 >= 0.0 && s3 >= 0.0) return true;
00231 return false;
00232 }
00233 }
00234 else if (dim_==3)
00235 {
00236 TEUCHOS_TEST_FOR_EXCEPT(true);
00237 return false;
00238 }
00239 else
00240 {
00241 TEUCHOS_TEST_FOR_EXCEPTION(dim_<=0 || dim_>3, std::runtime_error,
00242 "invalid point dimension " << dim_);
00243 return false;
00244 }
00245 }
00246
00247 bool AToCPointLocator::cellContainsPoint(int cellLID,
00248 const double* x,
00249 const int* facetLID,
00250 double* localCoords) const
00251 {
00252 if (dim_==2)
00253 {
00254 const double* A = mesh_.nodePositionView(facetLID[0]);
00255 const double* B = mesh_.nodePositionView(facetLID[1]);
00256 const double* C = mesh_.nodePositionView(facetLID[2]);
00257
00258
00259 double sign = orient2D(A, B, C);
00260 if (sign > 0.0)
00261 {
00262 double s1 = orient2D(A, B, x);
00263 double s2 = orient2D(B, C, x);
00264 double s3 = orient2D(C, A, x);
00265 if (s1 >= 0.0 && s2 >= 0.0 && s3 >= 0.0)
00266 {
00267 double bax = B[0] - A[0];
00268 double bay = B[1] - A[1];
00269 double cax = C[0] - A[0];
00270 double cay = C[1] - A[1];
00271 double delta = bax*cay - bay*cax;
00272
00273 double xax = x[0] - A[0];
00274 double xay = x[1] - A[1];
00275
00276 localCoords[0] = (cay*xax - cax*xay)/delta;
00277 localCoords[1] = (bax*xay - bay*xax)/delta;
00278 return true;
00279 }
00280 return false;
00281 }
00282 else
00283 {
00284 double s1 = orient2D(A, C, x);
00285 double s2 = orient2D(C, B, x);
00286 double s3 = orient2D(B, A, x);
00287 if (s1 >= 0.0 && s2 >= 0.0 && s3 >= 0.0)
00288 {
00289
00290 std::cout << "swapping!" << std::endl;
00291 double bax = C[0] - A[0];
00292 double bay = C[1] - A[1];
00293 double cax = B[0] - A[0];
00294 double cay = B[1] - A[1];
00295 double delta = bax*cay - bay*cax;
00296
00297 double xax = x[0] - A[0];
00298 double xay = x[1] - A[1];
00299
00300 localCoords[0] = (cay*xax - cax*xay)/delta;
00301 localCoords[1] = (bax*xay - bay*xax)/delta;
00302 return true;
00303 }
00304 return false;
00305 }
00306 }
00307 else if (dim_==3)
00308 {
00309 TEUCHOS_TEST_FOR_EXCEPT(true);
00310 return false;
00311 }
00312 else
00313 {
00314 TEUCHOS_TEST_FOR_EXCEPTION(dim_<=0 || dim_>3, std::runtime_error,
00315 "invalid point dimension " << dim_);
00316 return false;
00317 }
00318 }
00319
00320 int AToCPointLocator::findEnclosingCell(int initialGuessLID,
00321 const double* x) const
00322 {
00323 std::queue<int> Q;
00324 Set<int> repeats;
00325 int d = mesh_.spatialDim();
00326 Array<int> facets(mesh_.numFacets(d, 0, 0));
00327 Array<int> facetOr(facets.size());
00328
00329 Q.push(initialGuessLID);
00330
00331 while (!Q.empty())
00332 {
00333 int next = Q.front();
00334 Q.pop();
00335 if (repeats.contains(next)) continue;
00336
00337
00338 mesh_.getFacetArray(d, next, 0, facets, facetOr);
00339 if (cellContainsPoint(next, x, &(facets[0]))) return next;
00340 repeats.put(next);
00341
00342 fillMaximalNeighbors(next, &(facets[0]));
00343
00344 for (Set<int>::const_iterator
00345 i=neighborSet_[next]->begin(); i!=neighborSet_[next]->end(); i++)
00346 {
00347 Q.push(*i);
00348 }
00349 }
00350 return -1;
00351 }
00352
00353
00354
00355
00356 int AToCPointLocator::findEnclosingCell(int initialGuessLID,
00357 const double* x,
00358 double* xLocal) const
00359 {
00360 std::queue<int> Q;
00361 Set<int> repeats;
00362 int d = mesh_.spatialDim();
00363 Array<int> facets(mesh_.numFacets(d, 0, 0));
00364 Array<int> facetOr(facets.size());
00365
00366 Q.push(initialGuessLID);
00367
00368 while (!Q.empty())
00369 {
00370 int next = Q.front();
00371 Q.pop();
00372 if (repeats.contains(next)) continue;
00373
00374
00375 mesh_.getFacetArray(d, next, 0, facets, facetOr);
00376 if (cellContainsPoint(next, x, &(facets[0]), xLocal)) return next;
00377 repeats.put(next);
00378
00379 fillMaximalNeighbors(next, &(facets[0]));
00380
00381 for (Set<int>::const_iterator
00382 i=neighborSet_[next]->begin(); i!=neighborSet_[next]->end(); i++)
00383 {
00384 Q.push(*i);
00385 }
00386 }
00387 return -1;
00388 }
00389
00390
00391 Point AToCPointLocator::makePoint(int dim, const double* x)
00392 {
00393 if (dim==1) return Point(x[0]);
00394 else if (dim==2) return Point(x[0], x[1]);
00395 else return Point(x[0], x[1], x[2]);
00396 }