|
26 | 26 | #include <cassert> |
27 | 27 | #include <stdarg.h> |
28 | 28 |
|
| 29 | +static double PI = 3.14159265359; |
29 | 30 |
|
30 | 31 | namespace ma_dbg { |
31 | 32 |
|
@@ -256,5 +257,133 @@ void createCavityMesh(ma::Adapt* a, |
256 | 257 | createCavityMesh(a, tetsArray, prefix); |
257 | 258 | } |
258 | 259 |
|
| 260 | +static apf::Vector3 getPointOnEllipsoid( |
| 261 | + apf::Vector3 center, |
| 262 | + apf::Vector3 abc, |
| 263 | + apf::Matrix3x3 rotation, |
| 264 | + double scaleFactor, |
| 265 | + double u, |
| 266 | + double v) |
| 267 | +{ |
| 268 | + apf::Vector3 result; |
| 269 | + result[0] = abc[0] * cos(u) * cos(v); |
| 270 | + result[1] = abc[1] * cos(u) * sin(v); |
| 271 | + result[2] = abc[2] * sin(u); |
| 272 | + |
| 273 | + result = result * scaleFactor; |
| 274 | + |
| 275 | + result = rotation * result + center; |
| 276 | + return result; |
| 277 | +} |
| 278 | + |
| 279 | +static void makeEllipsoid( |
| 280 | + apf::Mesh2* msf, |
| 281 | + apf::Mesh2* mesh, |
| 282 | + apf::Field* sizes, |
| 283 | + apf::Field* frames, |
| 284 | + apf::MeshEntity* ent, |
| 285 | + int node, |
| 286 | + double scaleFactor, |
| 287 | + int sampleSize[2]) |
| 288 | +{ |
| 289 | + // first get the coordinate at node location |
| 290 | + apf::Vector3 xi; |
| 291 | + apf::Vector3 center; |
| 292 | + apf::FieldShape* fs = apf::getShape(sizes); |
| 293 | + fs->getNodeXi(mesh->getType(ent), node, xi); |
| 294 | + apf::MeshElement* me = apf::createMeshElement(mesh, ent); |
| 295 | + apf::mapLocalToGlobal(me, xi, center); |
| 296 | + apf::destroyMeshElement(me); |
| 297 | + |
| 298 | + |
| 299 | + // second get the sizes and frames at node |
| 300 | + apf::Vector3 abc; |
| 301 | + apf::getVector(sizes, ent, node, abc); |
| 302 | + |
| 303 | + apf::Matrix3x3 rotation; |
| 304 | + apf::getMatrix(frames, ent, node, rotation); |
| 305 | + |
| 306 | + |
| 307 | + double U0 = 0.0; |
| 308 | + double U1 = 2 * PI; |
| 309 | + double V0 = -PI/2.; |
| 310 | + double V1 = PI/2.; |
| 311 | + int n = sampleSize[0]; |
| 312 | + int m = sampleSize[1]; |
| 313 | + double dU = (U1 - U0) / (n-1); |
| 314 | + double dV = (V1 - V0) / (m-1); |
| 315 | + |
| 316 | + // make the array of vertex coordinates in the physical space |
| 317 | + std::vector<ma::Vector> ps; |
| 318 | + for (int j = 0; j < m; j++) { |
| 319 | + for (int i = 0; i < n; i++) { |
| 320 | + double u = U0 + i * dU; |
| 321 | + double v = V0 + j * dV; |
| 322 | + apf::Vector3 pt = getPointOnEllipsoid(center, abc, rotation, scaleFactor, u, v); |
| 323 | + ps.push_back(pt); |
| 324 | + } |
| 325 | + } |
| 326 | + // make the vertexes and set the coordinates using the array |
| 327 | + std::vector<apf::MeshEntity*> vs; |
| 328 | + for (size_t i = 0; i < ps.size(); i++) { |
| 329 | + apf::MeshEntity* newVert = msf->createVert(0); |
| 330 | + msf->setPoint(newVert, 0, ps[i]); |
| 331 | + vs.push_back(newVert); |
| 332 | + } |
| 333 | + |
| 334 | + PCU_ALWAYS_ASSERT(vs.size() == ps.size()); |
| 335 | + |
| 336 | + apf::MeshEntity* v[3]; |
| 337 | + // make the lower/upper t elems |
| 338 | + for (int i = 0; i < n-1; i++) { |
| 339 | + for (int j = 0; j < m-1; j++) { |
| 340 | + // upper triangle |
| 341 | + v[0] = vs[(i + 0) + n * (j + 0)]; |
| 342 | + v[1] = vs[(i + 0) + n * (j + 1)]; |
| 343 | + v[2] = vs[(i + 1) + n * (j + 0)]; |
| 344 | + apf::buildElement(msf, 0, apf::Mesh::TRIANGLE, v); |
| 345 | + // upper triangle |
| 346 | + v[0] = vs[(i + 0) + n * (j + 1)]; |
| 347 | + v[1] = vs[(i + 1) + n * (j + 1)]; |
| 348 | + v[2] = vs[(i + 1) + n * (j + 0)]; |
| 349 | + apf::buildElement(msf, 0, apf::Mesh::TRIANGLE, v); |
| 350 | + } |
| 351 | + } |
| 352 | +} |
| 353 | + |
| 354 | +void visualizeSizeField( |
| 355 | + apf::Mesh2* m, |
| 356 | + apf::Field* sizes, |
| 357 | + apf::Field* frames, |
| 358 | + int sampleSize[2], |
| 359 | + double userScale, |
| 360 | + const char* outputPrefix) |
| 361 | +{ |
| 362 | + // create the size-field visualization mesh |
| 363 | + apf::Mesh2* msf = apf::makeEmptyMdsMesh(gmi_load(".null"), 2, false); |
| 364 | + |
| 365 | + apf::FieldShape* fs = apf::getShape(sizes); |
| 366 | + int dim = m->getDimension(); |
| 367 | + |
| 368 | + apf::MeshEntity* ent; |
| 369 | + apf::MeshIterator* it; |
| 370 | + |
| 371 | + for (int d = 0; d <= dim; d++) { |
| 372 | + if (!fs->hasNodesIn(d)) continue; |
| 373 | + it = m->begin(d); |
| 374 | + while ( (ent = m->iterate(it)) ) { |
| 375 | + int type = m->getType(ent); |
| 376 | + int non = fs->countNodesOn(type); |
| 377 | + for (int n = 0; n < non; n++) |
| 378 | + makeEllipsoid(msf, m, sizes, frames, ent, n, userScale , sampleSize); |
| 379 | + } |
| 380 | + m->end(it); |
| 381 | + } |
| 382 | + |
| 383 | + apf::writeVtkFiles(outputPrefix, msf); |
| 384 | + |
| 385 | + msf->destroyNative(); |
| 386 | + apf::destroyMesh(msf); |
| 387 | +} |
259 | 388 |
|
260 | 389 | } |
0 commit comments