[370] | 1 | /* |
---|
| 2 | * ---------------------------------------------------------------------- |
---|
| 3 | * Rappture::MeshTri2D |
---|
| 4 | * This is a non-uniform, triangular mesh for 2-dimensional |
---|
| 5 | * structures. |
---|
| 6 | * |
---|
| 7 | * ====================================================================== |
---|
| 8 | * AUTHOR: Michael McLennan, Purdue University |
---|
| 9 | * Copyright (c) 2004-2006 Purdue Research Foundation |
---|
| 10 | * |
---|
| 11 | * See the file "license.terms" for information on usage and |
---|
| 12 | * redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES. |
---|
| 13 | * ====================================================================== |
---|
| 14 | */ |
---|
| 15 | #include <math.h> |
---|
[372] | 16 | #include <iostream> |
---|
[370] | 17 | #include "RpMeshTri2D.h" |
---|
| 18 | |
---|
| 19 | using namespace Rappture; |
---|
| 20 | |
---|
| 21 | CellTri2D::CellTri2D() |
---|
| 22 | { |
---|
| 23 | _cellId = -1; |
---|
| 24 | _nodes[0] = NULL; |
---|
| 25 | _nodes[1] = NULL; |
---|
| 26 | _nodes[2] = NULL; |
---|
| 27 | } |
---|
| 28 | |
---|
| 29 | CellTri2D::CellTri2D(int cellId, Node2D* n1Ptr, Node2D* n2Ptr, Node2D* n3Ptr) |
---|
| 30 | { |
---|
| 31 | _cellId = cellId; |
---|
| 32 | _nodes[0] = n1Ptr; |
---|
| 33 | _nodes[1] = n2Ptr; |
---|
| 34 | _nodes[2] = n3Ptr; |
---|
| 35 | } |
---|
| 36 | |
---|
| 37 | CellTri2D::CellTri2D(const CellTri2D& cell) |
---|
| 38 | { |
---|
| 39 | _cellId = cell._cellId; |
---|
| 40 | _nodes[0] = cell._nodes[0]; |
---|
| 41 | _nodes[1] = cell._nodes[1]; |
---|
| 42 | _nodes[2] = cell._nodes[2]; |
---|
| 43 | } |
---|
| 44 | |
---|
| 45 | CellTri2D& |
---|
| 46 | CellTri2D::operator=(const CellTri2D& cell) |
---|
| 47 | { |
---|
| 48 | _cellId = cell._cellId; |
---|
| 49 | _nodes[0] = cell._nodes[0]; |
---|
| 50 | _nodes[1] = cell._nodes[1]; |
---|
| 51 | _nodes[2] = cell._nodes[2]; |
---|
| 52 | return *this; |
---|
| 53 | } |
---|
| 54 | |
---|
| 55 | int |
---|
| 56 | CellTri2D::isNull() const |
---|
| 57 | { |
---|
| 58 | return (_nodes[0] == NULL || _nodes[1] == NULL || _nodes[2] == NULL); |
---|
| 59 | } |
---|
| 60 | |
---|
| 61 | int |
---|
| 62 | CellTri2D::isOutside() const |
---|
| 63 | { |
---|
| 64 | return (_cellId < 0); |
---|
| 65 | } |
---|
| 66 | |
---|
| 67 | void |
---|
| 68 | CellTri2D::clear() |
---|
| 69 | { |
---|
| 70 | _cellId = -1; |
---|
| 71 | _nodes[0] = _nodes[1] = _nodes[2] = NULL; |
---|
| 72 | } |
---|
| 73 | |
---|
| 74 | int |
---|
| 75 | CellTri2D::cellId() const |
---|
| 76 | { |
---|
| 77 | return _cellId; |
---|
| 78 | } |
---|
| 79 | |
---|
| 80 | int |
---|
| 81 | CellTri2D::nodeId(int n) const |
---|
| 82 | { |
---|
| 83 | assert(n >= 0 && n <= 2); |
---|
| 84 | return _nodes[n]->id(); |
---|
| 85 | } |
---|
| 86 | |
---|
| 87 | double |
---|
| 88 | CellTri2D::x(int n) const |
---|
| 89 | { |
---|
| 90 | assert(n >= 0 && n <= 2); |
---|
| 91 | return _nodes[n]->x(); |
---|
| 92 | } |
---|
| 93 | |
---|
| 94 | double |
---|
| 95 | CellTri2D::y(int n) const |
---|
| 96 | { |
---|
| 97 | assert(n >= 0 && n <= 2); |
---|
| 98 | return _nodes[n]->y(); |
---|
| 99 | } |
---|
| 100 | |
---|
| 101 | void |
---|
| 102 | CellTri2D::barycentrics(const Node2D& node, double* phi) const |
---|
| 103 | { |
---|
| 104 | assert( _nodes[0] != NULL && _nodes[1] != NULL && _nodes[2] != NULL); |
---|
| 105 | |
---|
| 106 | double x2 = _nodes[1]->x() - _nodes[0]->x(); |
---|
| 107 | double y2 = _nodes[1]->y() - _nodes[0]->y(); |
---|
| 108 | double x3 = _nodes[2]->x() - _nodes[0]->x(); |
---|
| 109 | double y3 = _nodes[2]->y() - _nodes[0]->y(); |
---|
| 110 | double xr = node.x() - _nodes[0]->x(); |
---|
| 111 | double yr = node.y() - _nodes[0]->y(); |
---|
| 112 | double det = x2*y3-x3*y2; |
---|
| 113 | |
---|
[372] | 114 | if (det != 0.0) { |
---|
| 115 | phi[1] = (xr*y3 - x3*yr)/det; |
---|
| 116 | phi[2] = (x2*yr - xr*y2)/det; |
---|
| 117 | phi[0] = 1.0-phi[1]-phi[2]; |
---|
| 118 | } else { |
---|
| 119 | phi[1] = NAN; |
---|
| 120 | phi[2] = NAN; |
---|
| 121 | phi[0] = NAN; |
---|
| 122 | } |
---|
[370] | 123 | } |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | Tri2D::Tri2D() |
---|
| 127 | { |
---|
| 128 | nodes[0] = nodes[1] = nodes[2] = -1; |
---|
| 129 | neighbors[0] = neighbors[1] = neighbors[2] = -1; |
---|
| 130 | } |
---|
| 131 | |
---|
| 132 | Tri2D::Tri2D(int n1, int n2, int n3) |
---|
| 133 | { |
---|
| 134 | nodes[0] = n1; |
---|
| 135 | nodes[1] = n2; |
---|
| 136 | nodes[2] = n3; |
---|
| 137 | neighbors[0] = neighbors[1] = neighbors[2] = -1; |
---|
| 138 | } |
---|
| 139 | |
---|
| 140 | |
---|
| 141 | MeshTri2D::MeshTri2D() |
---|
| 142 | : _counter(0), |
---|
| 143 | _id2nodeDirty(0), |
---|
| 144 | _id2node(100,-1) |
---|
| 145 | { |
---|
| 146 | _nodelist.reserve(1024); |
---|
| 147 | _min[0] = _min[1] = NAN; |
---|
| 148 | _max[0] = _max[1] = NAN; |
---|
| 149 | } |
---|
| 150 | |
---|
| 151 | MeshTri2D::MeshTri2D(const MeshTri2D& mesh) |
---|
| 152 | : _nodelist(mesh._nodelist), |
---|
| 153 | _counter(mesh._counter), |
---|
| 154 | _celllist(mesh._celllist), |
---|
| 155 | _id2nodeDirty(mesh._id2nodeDirty), |
---|
| 156 | _id2node(mesh._id2node) |
---|
| 157 | { |
---|
| 158 | for (int i=0; i < 2; i++) { |
---|
| 159 | _min[i] = mesh._min[i]; |
---|
| 160 | _max[i] = mesh._max[i]; |
---|
| 161 | } |
---|
| 162 | } |
---|
| 163 | |
---|
| 164 | MeshTri2D& |
---|
| 165 | MeshTri2D::operator=(const MeshTri2D& mesh) |
---|
| 166 | { |
---|
| 167 | _nodelist = mesh._nodelist; |
---|
| 168 | _counter = mesh._counter; |
---|
| 169 | for (int i=0; i < 2; i++) { |
---|
| 170 | _min[i] = mesh._min[i]; |
---|
| 171 | _max[i] = mesh._max[i]; |
---|
| 172 | } |
---|
| 173 | _celllist = mesh._celllist; |
---|
| 174 | _id2nodeDirty = mesh._id2nodeDirty; |
---|
| 175 | _id2node = mesh._id2node; |
---|
| 176 | _lastLocate.clear(); |
---|
| 177 | } |
---|
| 178 | |
---|
| 179 | MeshTri2D::~MeshTri2D() |
---|
| 180 | { |
---|
| 181 | } |
---|
| 182 | |
---|
| 183 | Node2D& |
---|
| 184 | MeshTri2D::addNode(const Node2D& nd) |
---|
| 185 | { |
---|
| 186 | Node2D node(nd); |
---|
| 187 | |
---|
| 188 | if (node.id() < 0) { |
---|
| 189 | node.id(_counter++); |
---|
| 190 | } else { |
---|
| 191 | // see if this node already exists |
---|
| 192 | Node2D* nptr = _getNodeById(node.id()); |
---|
| 193 | if (nptr) { |
---|
| 194 | return *nptr; |
---|
| 195 | } |
---|
| 196 | } |
---|
| 197 | |
---|
| 198 | // add this new node |
---|
| 199 | _nodelist.push_back(node); |
---|
| 200 | |
---|
| 201 | if (isnan(_min[0])) { |
---|
| 202 | _min[0] = _max[0] = node.x(); |
---|
| 203 | _min[1] = _max[1] = node.y(); |
---|
| 204 | } else { |
---|
| 205 | if (node.x() < _min[0]) { _min[0] = node.x(); } |
---|
| 206 | if (node.x() > _max[0]) { _max[0] = node.x(); } |
---|
| 207 | if (node.y() < _min[1]) { _min[1] = node.y(); } |
---|
| 208 | if (node.y() > _max[1]) { _max[1] = node.y(); } |
---|
| 209 | } |
---|
| 210 | |
---|
| 211 | // index of this new node |
---|
| 212 | int n = _nodelist.size()-1; |
---|
| 213 | |
---|
| 214 | if (!_id2nodeDirty) { |
---|
| 215 | // id2node map up to date? then keep it up to date |
---|
| 216 | if (node.id() >= _id2node.size()) { |
---|
| 217 | int newsize = 2*_id2node.size(); |
---|
| 218 | _id2node.reserve(newsize); |
---|
| 219 | for (int i=_id2node.size(); i < newsize; i++) { |
---|
| 220 | _id2node[i] = -1; |
---|
| 221 | } |
---|
| 222 | } |
---|
| 223 | _id2node[node.id()] = n; |
---|
| 224 | } |
---|
| 225 | return _nodelist[n]; |
---|
| 226 | } |
---|
| 227 | |
---|
| 228 | void |
---|
| 229 | MeshTri2D::addCell(const Node2D& n1, const Node2D& n2, const Node2D& n3) |
---|
| 230 | { |
---|
| 231 | Node2D node1 = addNode(n1); |
---|
| 232 | Node2D node2 = addNode(n2); |
---|
| 233 | Node2D node3 = addNode(n3); |
---|
| 234 | addCell(node1.id(), node2.id(), node3.id()); |
---|
| 235 | } |
---|
| 236 | |
---|
| 237 | void |
---|
| 238 | MeshTri2D::addCell(int nId1, int nId2, int nId3) |
---|
| 239 | { |
---|
| 240 | _celllist.push_back( Tri2D(nId1,nId2,nId3) ); |
---|
| 241 | int triId = _celllist.size()-1; |
---|
| 242 | |
---|
| 243 | Edge2D edge; |
---|
| 244 | int nodes[4]; |
---|
| 245 | nodes[0] = nId1; |
---|
| 246 | nodes[1] = nId2; |
---|
| 247 | nodes[2] = nId3; |
---|
| 248 | nodes[3] = nId1; |
---|
| 249 | |
---|
| 250 | // update the neighbors for this triangle |
---|
| 251 | for (int i=0; i < 3; i++) { |
---|
| 252 | int n = (i+2) % 3; // this is node opposite i/i+1 edge |
---|
| 253 | |
---|
| 254 | // build an edge with nodes in proper order |
---|
| 255 | if (nodes[i] < nodes[i+1]) { |
---|
| 256 | edge.fromNode = nodes[i]; |
---|
| 257 | edge.toNode = nodes[i+1]; |
---|
| 258 | } else { |
---|
| 259 | edge.fromNode = nodes[i+1]; |
---|
| 260 | edge.toNode = nodes[i]; |
---|
| 261 | } |
---|
| 262 | |
---|
| 263 | Neighbor2D& nbr = _edge2neighbor[edge]; |
---|
| 264 | if (nbr.triId < 0) { |
---|
| 265 | // not found? then register this for later |
---|
| 266 | nbr.triId = triId; |
---|
| 267 | nbr.index = n; |
---|
| 268 | } else { |
---|
| 269 | // this triangle points to other one for the edge |
---|
| 270 | _celllist[triId].neighbors[n] = nbr.triId; |
---|
| 271 | // other triangle points to this one for the same edge |
---|
| 272 | _celllist[nbr.triId].neighbors[nbr.index] = triId; |
---|
| 273 | _edge2neighbor.erase(edge); |
---|
| 274 | } |
---|
| 275 | } |
---|
| 276 | } |
---|
| 277 | |
---|
| 278 | MeshTri2D& |
---|
| 279 | MeshTri2D::clear() |
---|
| 280 | { |
---|
| 281 | _nodelist.clear(); |
---|
| 282 | _counter = 0; |
---|
| 283 | _id2nodeDirty = 0; |
---|
| 284 | _id2node.assign(100, -1); |
---|
| 285 | _lastLocate.clear(); |
---|
| 286 | } |
---|
| 287 | |
---|
| 288 | int |
---|
| 289 | MeshTri2D::sizeNodes() const |
---|
| 290 | { |
---|
| 291 | return _nodelist.size(); |
---|
| 292 | } |
---|
| 293 | |
---|
| 294 | Node2D& |
---|
| 295 | MeshTri2D::atNode(int pos) |
---|
| 296 | { |
---|
| 297 | assert(pos >= 0 && pos < _nodelist.size()); |
---|
| 298 | return _nodelist.at(pos); |
---|
| 299 | } |
---|
| 300 | |
---|
| 301 | int |
---|
| 302 | MeshTri2D::sizeCells() const |
---|
| 303 | { |
---|
| 304 | return _celllist.size(); |
---|
| 305 | } |
---|
| 306 | |
---|
| 307 | CellTri2D |
---|
| 308 | MeshTri2D::atCell(int pos) |
---|
| 309 | { |
---|
| 310 | assert(pos >= 0 && pos < _celllist.size()); |
---|
| 311 | |
---|
| 312 | Tri2D& cell = _celllist[pos]; |
---|
| 313 | Node2D* n1Ptr = _getNodeById(cell.nodes[0]); |
---|
| 314 | Node2D* n2Ptr = _getNodeById(cell.nodes[1]); |
---|
| 315 | Node2D* n3Ptr = _getNodeById(cell.nodes[2]); |
---|
| 316 | assert(n1Ptr && n2Ptr && n3Ptr); |
---|
| 317 | |
---|
| 318 | CellTri2D rval(pos, n1Ptr, n2Ptr, n3Ptr); |
---|
| 319 | return rval; |
---|
| 320 | } |
---|
| 321 | |
---|
| 322 | double |
---|
| 323 | MeshTri2D::rangeMin(Axis which) const |
---|
| 324 | { |
---|
| 325 | assert(which != Rappture::zaxis); |
---|
| 326 | return _min[which]; |
---|
| 327 | } |
---|
| 328 | |
---|
| 329 | double |
---|
| 330 | MeshTri2D::rangeMax(Axis which) const |
---|
| 331 | { |
---|
| 332 | assert(which != Rappture::zaxis); |
---|
| 333 | return _max[which]; |
---|
| 334 | } |
---|
| 335 | |
---|
| 336 | CellTri2D |
---|
| 337 | MeshTri2D::locate(const Node2D& node) const |
---|
| 338 | { |
---|
| 339 | MeshTri2D* nonconst = (MeshTri2D*)this; |
---|
| 340 | CellTri2D cell = _lastLocate; |
---|
| 341 | |
---|
| 342 | if (cell.isNull() && _celllist.size() > 0) { |
---|
| 343 | Tri2D& tri = nonconst->_celllist[0]; |
---|
| 344 | cell = CellTri2D(0, &nonconst->_nodelist[tri.nodes[0]], |
---|
| 345 | &nonconst->_nodelist[tri.nodes[1]], |
---|
| 346 | &nonconst->_nodelist[tri.nodes[2]]); |
---|
| 347 | } |
---|
| 348 | |
---|
[372] | 349 | int lastCellId = -1; |
---|
[370] | 350 | while (!cell.isNull()) { |
---|
| 351 | double phi[3]; |
---|
| 352 | |
---|
| 353 | // compute barycentric coords |
---|
| 354 | // if all are >= 0, then this tri contains node |
---|
| 355 | cell.barycentrics(node, phi); |
---|
[372] | 356 | if (isnan(phi[0])) { |
---|
| 357 | cell.clear(); |
---|
| 358 | return cell; |
---|
| 359 | } |
---|
[370] | 360 | if (phi[0] >= 0.0 && phi[1] >= 0.0 && phi[2] >= 0.0) { |
---|
| 361 | break; |
---|
| 362 | } |
---|
| 363 | |
---|
| 364 | // find the smallest (most negative) coord phi, and search that dir |
---|
| 365 | int dir = 0; |
---|
| 366 | for (int i=1; i <= 2; i++) { |
---|
| 367 | if (phi[i] < phi[dir]) { |
---|
| 368 | dir = i; |
---|
| 369 | } |
---|
| 370 | } |
---|
| 371 | |
---|
| 372 | Tri2D& tri = nonconst->_celllist[ cell.cellId() ]; |
---|
| 373 | int neighborId = tri.neighbors[dir]; |
---|
[372] | 374 | if (neighborId < 0 || neighborId == lastCellId) { |
---|
| 375 | if (neighborId == lastCellId) { |
---|
| 376 | std::cout << "loop at " << lastCellId << "->" << neighborId << std::endl; |
---|
| 377 | } |
---|
[370] | 378 | cell.clear(); |
---|
| 379 | return cell; |
---|
| 380 | } |
---|
| 381 | |
---|
[372] | 382 | lastCellId = cell.cellId(); |
---|
[370] | 383 | Tri2D& tri2 = nonconst->_celllist[neighborId]; |
---|
| 384 | Node2D *n1Ptr = &nonconst->_nodelist[tri2.nodes[0]]; |
---|
| 385 | Node2D *n2Ptr = &nonconst->_nodelist[tri2.nodes[1]]; |
---|
| 386 | Node2D *n3Ptr = &nonconst->_nodelist[tri2.nodes[2]]; |
---|
| 387 | cell = CellTri2D(neighborId, n1Ptr, n2Ptr, n3Ptr); |
---|
| 388 | } |
---|
| 389 | |
---|
| 390 | nonconst->_lastLocate = cell; |
---|
| 391 | return cell; |
---|
| 392 | } |
---|
| 393 | |
---|
| 394 | Ptr<Serializable> |
---|
| 395 | MeshTri2D::create() |
---|
| 396 | { |
---|
| 397 | return Ptr<Serializable>( (Serializable*) new MeshTri2D() ); |
---|
| 398 | } |
---|
| 399 | |
---|
| 400 | void |
---|
| 401 | MeshTri2D::serialize_A(SerialBuffer& buffer) const |
---|
| 402 | { |
---|
| 403 | } |
---|
| 404 | |
---|
| 405 | Outcome |
---|
| 406 | deserialize_A(SerialBuffer& buffer) |
---|
| 407 | { |
---|
| 408 | Outcome status; |
---|
| 409 | return status; |
---|
| 410 | } |
---|
| 411 | |
---|
| 412 | Node2D* |
---|
| 413 | MeshTri2D::_getNodeById(int nodeId) |
---|
| 414 | { |
---|
| 415 | MeshTri2D* nonconst = (MeshTri2D*)this; |
---|
| 416 | Node2D *rptr = NULL; |
---|
| 417 | nonconst->_rebuildNodeIdMap(); |
---|
| 418 | |
---|
| 419 | if (nodeId < _id2node.size()) { |
---|
| 420 | int n = _id2node[nodeId]; |
---|
| 421 | if (n >= 0) { |
---|
| 422 | return &_nodelist[n]; |
---|
| 423 | } |
---|
| 424 | } |
---|
| 425 | return rptr; |
---|
| 426 | } |
---|
| 427 | |
---|
| 428 | void |
---|
| 429 | MeshTri2D::_rebuildNodeIdMap() |
---|
| 430 | { |
---|
| 431 | if (_id2nodeDirty) { |
---|
| 432 | _min[0] = _min[1] = NAN; |
---|
| 433 | _max[0] = _max[1] = NAN; |
---|
| 434 | |
---|
| 435 | // figure out how big the _id2node array should be |
---|
| 436 | int maxId = -1; |
---|
| 437 | std::vector<Node2D>::iterator n = _nodelist.begin(); |
---|
| 438 | while (n != _nodelist.end()) { |
---|
| 439 | if (n->id() > maxId) { |
---|
| 440 | maxId = n->id(); |
---|
| 441 | } |
---|
| 442 | ++n; |
---|
| 443 | } |
---|
| 444 | if (maxId > 0) { |
---|
| 445 | _id2node.assign(maxId+1, -1); |
---|
| 446 | } |
---|
| 447 | |
---|
| 448 | // scan through and map id -> node index |
---|
| 449 | int i = 0; |
---|
| 450 | n = _nodelist.begin(); |
---|
| 451 | while (n != _nodelist.end()) { |
---|
| 452 | _id2node[n->id()] = i++; |
---|
| 453 | |
---|
| 454 | if (isnan(_min[0])) { |
---|
| 455 | _min[0] = _max[0] = n->x(); |
---|
| 456 | _min[1] = _max[1] = n->y(); |
---|
| 457 | } else { |
---|
| 458 | if (n->x() < _min[0]) { _min[0] = n->x(); } |
---|
| 459 | if (n->x() > _max[0]) { _max[0] = n->x(); } |
---|
| 460 | if (n->y() < _min[1]) { _min[1] = n->y(); } |
---|
| 461 | if (n->y() > _max[1]) { _max[1] = n->y(); } |
---|
| 462 | } |
---|
| 463 | ++n; |
---|
| 464 | } |
---|
| 465 | _id2nodeDirty = 0; |
---|
| 466 | } |
---|
| 467 | } |
---|