Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat (limited to 'extern/mantaflow/preprocessed/plugin/meshplugins.cpp')
-rw-r--r--extern/mantaflow/preprocessed/plugin/meshplugins.cpp780
1 files changed, 780 insertions, 0 deletions
diff --git a/extern/mantaflow/preprocessed/plugin/meshplugins.cpp b/extern/mantaflow/preprocessed/plugin/meshplugins.cpp
new file mode 100644
index 00000000000..415bca153d0
--- /dev/null
+++ b/extern/mantaflow/preprocessed/plugin/meshplugins.cpp
@@ -0,0 +1,780 @@
+
+
+// DO NOT EDIT !
+// This file is generated using the MantaFlow preprocessor (prep generate).
+
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2011 Tobias Pfaff, Nils Thuerey
+ *
+ * This program is free software, distributed under the terms of the
+ * Apache License, Version 2.0
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Smoothing etc. for meshes
+ *
+ ******************************************************************************/
+
+/******************************************************************************/
+// Copyright note:
+//
+// These functions (C) Chris Wojtan
+// Long-term goal is to unify with his split&merge codebase
+//
+/******************************************************************************/
+
+#include <queue>
+#include <algorithm>
+#include "mesh.h"
+#include "kernel.h"
+#include "edgecollapse.h"
+#include <mesh.h>
+#include <stack>
+
+using namespace std;
+
+namespace Manta {
+
+//! Mesh smoothing
+/*! see Desbrun 99 "Implicit fairing of of irregular meshes using diffusion and curvature flow"*/
+void smoothMesh(Mesh &mesh, Real strength, int steps = 1, Real minLength = 1e-5)
+{
+ const Real dt = mesh.getParent()->getDt();
+ const Real str = min(dt * strength, (Real)1);
+ mesh.rebuildQuickCheck();
+
+ // calculate original mesh volume
+ Vec3 origCM;
+ Real origVolume = mesh.computeCenterOfMass(origCM);
+
+ // temp vertices
+ const int numCorners = mesh.numTris() * 3;
+ const int numNodes = mesh.numNodes();
+ vector<Vec3> temp(numNodes);
+ vector<bool> visited(numNodes);
+
+ for (int s = 0; s < steps; s++) {
+ // reset markers
+ for (size_t i = 0; i < visited.size(); i++)
+ visited[i] = false;
+
+ for (int c = 0; c < numCorners; c++) {
+ const int node = mesh.corners(c).node;
+ if (visited[node])
+ continue;
+
+ const Vec3 pos = mesh.nodes(node).pos;
+ Vec3 dx(0.0);
+ Real totalLen = 0;
+
+ // rotate around vertex
+ set<int> &ring = mesh.get1Ring(node).nodes;
+ for (set<int>::iterator it = ring.begin(); it != ring.end(); it++) {
+ Vec3 edge = mesh.nodes(*it).pos - pos;
+ Real len = norm(edge);
+
+ if (len > minLength) {
+ dx += edge * (1.0 / len);
+ totalLen += len;
+ }
+ else {
+ totalLen = 0.0;
+ break;
+ }
+ }
+ visited[node] = true;
+ temp[node] = pos;
+ if (totalLen != 0)
+ temp[node] += dx * (str / totalLen);
+ }
+
+ // copy back
+ for (int n = 0; n < numNodes; n++)
+ if (!mesh.isNodeFixed(n))
+ mesh.nodes(n).pos = temp[n];
+ }
+
+ // calculate new mesh volume
+ Vec3 newCM;
+ Real newVolume = mesh.computeCenterOfMass(newCM);
+
+ // preserve volume : scale relative to CM
+ Real beta;
+#if defined(WIN32) || defined(_WIN32)
+ beta = pow((Real)std::abs(origVolume / newVolume), (Real)(1. / 3.));
+#else
+ beta = cbrt(origVolume / newVolume);
+#endif
+
+ for (int n = 0; n < numNodes; n++)
+ if (!mesh.isNodeFixed(n))
+ mesh.nodes(n).pos = origCM + (mesh.nodes(n).pos - newCM) * beta;
+}
+static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+{
+ try {
+ PbArgs _args(_linargs, _kwds);
+ FluidSolver *parent = _args.obtainParent();
+ bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+ pbPreparePlugin(parent, "smoothMesh", !noTiming);
+ PyObject *_retval = 0;
+ {
+ ArgLocker _lock;
+ Mesh &mesh = *_args.getPtr<Mesh>("mesh", 0, &_lock);
+ Real strength = _args.get<Real>("strength", 1, &_lock);
+ int steps = _args.getOpt<int>("steps", 2, 1, &_lock);
+ Real minLength = _args.getOpt<Real>("minLength", 3, 1e-5, &_lock);
+ _retval = getPyNone();
+ smoothMesh(mesh, strength, steps, minLength);
+ _args.check();
+ }
+ pbFinalizePlugin(parent, "smoothMesh", !noTiming);
+ return _retval;
+ }
+ catch (std::exception &e) {
+ pbSetError("smoothMesh", e.what());
+ return 0;
+ }
+}
+static const Pb::Register _RP_smoothMesh("", "smoothMesh", _W_0);
+extern "C" {
+void PbRegister_smoothMesh()
+{
+ KEEP_UNUSED(_RP_smoothMesh);
+}
+}
+
+//! Subdivide and edgecollapse to guarantee mesh with edgelengths between
+//! min/maxLength and an angle below minAngle
+void subdivideMesh(
+ Mesh &mesh, Real minAngle, Real minLength, Real maxLength, bool cutTubes = false)
+{
+ // gather some statistics
+ int edgeSubdivs = 0, edgeCollsAngle = 0, edgeCollsLen = 0, edgeKill = 0;
+ mesh.rebuildQuickCheck();
+
+ vector<int> deletedNodes;
+ map<int, bool> taintedTris;
+ priority_queue<pair<Real, int>> pq;
+
+ //////////////////////////////////////////
+ // EDGE COLLAPSE //
+ // - particles marked for deletation //
+ //////////////////////////////////////////
+
+ for (int t = 0; t < mesh.numTris(); t++) {
+ if (taintedTris.find(t) != taintedTris.end())
+ continue;
+
+ // check if at least 2 nodes are marked for delete
+ bool k[3];
+ int numKill = 0;
+ for (int i = 0; i < 3; i++) {
+ k[i] = mesh.nodes(mesh.tris(t).c[i]).flags & Mesh::NfKillme;
+ if (k[i])
+ numKill++;
+ }
+ if (numKill < 2)
+ continue;
+
+ if (k[0] && k[1])
+ CollapseEdge(mesh,
+ t,
+ 2,
+ mesh.getEdge(t, 0),
+ mesh.getNode(t, 0),
+ deletedNodes,
+ taintedTris,
+ edgeKill,
+ cutTubes);
+ else if (k[1] && k[2])
+ CollapseEdge(mesh,
+ t,
+ 0,
+ mesh.getEdge(t, 1),
+ mesh.getNode(t, 1),
+ deletedNodes,
+ taintedTris,
+ edgeKill,
+ cutTubes);
+ else if (k[2] && k[0])
+ CollapseEdge(mesh,
+ t,
+ 1,
+ mesh.getEdge(t, 2),
+ mesh.getNode(t, 2),
+ deletedNodes,
+ taintedTris,
+ edgeKill,
+ cutTubes);
+ }
+
+ //////////////////////////////////////////
+ // EDGE COLLAPSING //
+ // - based on small triangle angle //
+ //////////////////////////////////////////
+
+ if (minAngle > 0) {
+ for (int t = 0; t < mesh.numTris(); t++) {
+ // we only want to run through the edge list ONCE.
+ // we achieve this in a method very similar to the above subdivision method.
+
+ // if this triangle has already been deleted, ignore it
+ if (taintedTris.find(t) != taintedTris.end())
+ continue;
+
+ // first we find the angles of this triangle
+ Vec3 e0 = mesh.getEdge(t, 0), e1 = mesh.getEdge(t, 1), e2 = mesh.getEdge(t, 2);
+ Vec3 ne0 = e0;
+ Vec3 ne1 = e1;
+ Vec3 ne2 = e2;
+ normalize(ne0);
+ normalize(ne1);
+ normalize(ne2);
+
+ // Real thisArea = sqrMag(cross(-e2,e0));
+ // small angle approximation says sin(x) = arcsin(x) = x,
+ // arccos(x) = pi/2 - arcsin(x),
+ // cos(x) = dot(A,B),
+ // so angle is approximately 1 - dot(A,B).
+ Real angle[3];
+ angle[0] = 1.0 - dot(ne0, -ne2);
+ angle[1] = 1.0 - dot(ne1, -ne0);
+ angle[2] = 1.0 - dot(ne2, -ne1);
+ Real worstAngle = angle[0];
+ int which = 0;
+ if (angle[1] < worstAngle) {
+ worstAngle = angle[1];
+ which = 1;
+ }
+ if (angle[2] < worstAngle) {
+ worstAngle = angle[2];
+ which = 2;
+ }
+
+ // then we see if the angle is too small
+ if (worstAngle < minAngle) {
+ Vec3 edgevect;
+ Vec3 endpoint;
+ switch (which) {
+ case 0:
+ endpoint = mesh.getNode(t, 1);
+ edgevect = e1;
+ break;
+ case 1:
+ endpoint = mesh.getNode(t, 2);
+ edgevect = e2;
+ break;
+ case 2:
+ endpoint = mesh.getNode(t, 0);
+ edgevect = e0;
+ break;
+ default:
+ break;
+ }
+
+ CollapseEdge(mesh,
+ t,
+ which,
+ edgevect,
+ endpoint,
+ deletedNodes,
+ taintedTris,
+ edgeCollsAngle,
+ cutTubes);
+ }
+ }
+ }
+
+ //////////////////////
+ // EDGE SUBDIVISION //
+ //////////////////////
+
+ Real maxLength2 = maxLength * maxLength;
+ for (int t = 0; t < mesh.numTris(); t++) {
+ // first we find the maximum length edge in this triangle
+ Vec3 e0 = mesh.getEdge(t, 0), e1 = mesh.getEdge(t, 1), e2 = mesh.getEdge(t, 2);
+ Real d0 = normSquare(e0);
+ Real d1 = normSquare(e1);
+ Real d2 = normSquare(e2);
+
+ Real longest = max(d0, max(d1, d2));
+ if (longest > maxLength2) {
+ pq.push(pair<Real, int>(longest, t));
+ }
+ }
+ if (maxLength > 0) {
+
+ while (!pq.empty() && pq.top().first > maxLength2) {
+ // we only want to run through the edge list ONCE
+ // and we want to subdivide the original edges before we subdivide any newer, shorter edges,
+ // so whenever we subdivide, we add the 2 new triangles on the end of the SurfaceTri vector
+ // and mark the original subdivided triangles for deletion.
+ // when we are done subdividing, we delete the obsolete triangles
+
+ int triA = pq.top().second;
+ pq.pop();
+
+ if (taintedTris.find(triA) != taintedTris.end())
+ continue;
+
+ // first we find the maximum length edge in this triangle
+ Vec3 e0 = mesh.getEdge(triA, 0), e1 = mesh.getEdge(triA, 1), e2 = mesh.getEdge(triA, 2);
+ Real d0 = normSquare(e0);
+ Real d1 = normSquare(e1);
+ Real d2 = normSquare(e2);
+
+ Vec3 edgevect;
+ Vec3 endpoint;
+ int which;
+ if (d0 > d1) {
+ if (d0 > d2) {
+ edgevect = e0;
+ endpoint = mesh.getNode(triA, 0);
+ ;
+ which = 2; // 2 opposite of edge 0-1
+ }
+ else {
+ edgevect = e2;
+ endpoint = mesh.getNode(triA, 2);
+ which = 1; // 1 opposite of edge 2-0
+ }
+ }
+ else {
+ if (d1 > d2) {
+ edgevect = e1;
+ endpoint = mesh.getNode(triA, 1);
+ which = 0; // 0 opposite of edge 1-2
+ }
+ else {
+ edgevect = e2;
+ endpoint = mesh.getNode(triA, 2);
+ which = 1; // 1 opposite of edge 2-0
+ }
+ }
+ // This edge is too long, so we split it in the middle
+
+ // *
+ // / \.
+ // /C0 \.
+ // / \.
+ // / \.
+ // / B \.
+ // / \.
+ // /C1 C2 \.
+ // *---------------*
+ // \C2 C1 /
+ // \ /
+ // \ A /
+ // \ /
+ // \ /
+ // \C0 /
+ // \ /
+ // *
+ //
+ // BECOMES
+ //
+ // *
+ // /|\.
+ // / | \.
+ // /C0|C0\.
+ // / | \.
+ // / B1 | B2 \.
+ // / | \.
+ // /C1 C2|C1 C2 \.
+ // *-------*-------*
+ // \C2 C1|C2 C1/
+ // \ | /
+ // \ A2 | A1 /
+ // \ | /
+ // \C0|C0/
+ // \ | /
+ // \|/
+ // *
+
+ int triB = -1;
+ bool haveB = false;
+ Corner ca_old[3], cb_old[3];
+ ca_old[0] = mesh.corners(triA, which);
+ ca_old[1] = mesh.corners(ca_old[0].next);
+ ca_old[2] = mesh.corners(ca_old[0].prev);
+ if (ca_old[0].opposite >= 0) {
+ cb_old[0] = mesh.corners(ca_old[0].opposite);
+ cb_old[1] = mesh.corners(cb_old[0].next);
+ cb_old[2] = mesh.corners(cb_old[0].prev);
+ triB = cb_old[0].tri;
+ haveB = true;
+ }
+ // else throw Error("nonmanifold");
+
+ // subdivide in the middle of the edge and create new triangles
+ Node newNode;
+ newNode.flags = 0;
+
+ newNode.pos = endpoint + 0.5 * edgevect; // fallback: linear average
+ // default: use butterfly
+ if (haveB)
+ newNode.pos = ModifiedButterflySubdivision(mesh, ca_old[0], cb_old[0], newNode.pos);
+
+ // find indices of two points of 'which'-edge
+ // merge flags
+ int P0 = ca_old[1].node;
+ int P1 = ca_old[2].node;
+ newNode.flags = mesh.nodes(P0).flags | mesh.nodes(P1).flags;
+
+ Real len0 = norm(mesh.nodes(P0).pos - newNode.pos);
+ Real len1 = norm(mesh.nodes(P1).pos - newNode.pos);
+
+ // remove P0/P1 1-ring connection
+ mesh.get1Ring(P0).nodes.erase(P1);
+ mesh.get1Ring(P1).nodes.erase(P0);
+ mesh.get1Ring(P0).tris.erase(triA);
+ mesh.get1Ring(P1).tris.erase(triA);
+ mesh.get1Ring(ca_old[0].node).tris.erase(triA);
+ if (haveB) {
+ mesh.get1Ring(P0).tris.erase(triB);
+ mesh.get1Ring(P1).tris.erase(triB);
+ mesh.get1Ring(cb_old[0].node).tris.erase(triB);
+ }
+
+ // init channel properties for new node
+ for (int i = 0; i < mesh.numNodeChannels(); i++) {
+ mesh.nodeChannel(i)->addInterpol(P0, P1, len0 / (len0 + len1));
+ }
+
+ // write to array
+ mesh.addTri(Triangle(ca_old[0].node, ca_old[1].node, mesh.numNodes()));
+ mesh.addTri(Triangle(ca_old[0].node, mesh.numNodes(), ca_old[2].node));
+ if (haveB) {
+ mesh.addTri(Triangle(cb_old[0].node, cb_old[1].node, mesh.numNodes()));
+ mesh.addTri(Triangle(cb_old[0].node, mesh.numNodes(), cb_old[2].node));
+ }
+ mesh.addNode(newNode);
+
+ const int nt = haveB ? 4 : 2;
+ int triA1 = mesh.numTris() - nt;
+ int triA2 = mesh.numTris() - nt + 1;
+ int triB1 = 0, triB2 = 0;
+ if (haveB) {
+ triB1 = mesh.numTris() - nt + 2;
+ triB2 = mesh.numTris() - nt + 3;
+ }
+ mesh.tris(triA1).flags = mesh.tris(triA).flags;
+ mesh.tris(triA2).flags = mesh.tris(triA).flags;
+ mesh.tris(triB1).flags = mesh.tris(triB).flags;
+ mesh.tris(triB2).flags = mesh.tris(triB).flags;
+
+ // connect new triangles to outside triangles,
+ // and connect outside triangles to these new ones
+ for (int c = 0; c < 3; c++)
+ mesh.addCorner(Corner(triA1, mesh.tris(triA1).c[c]));
+ for (int c = 0; c < 3; c++)
+ mesh.addCorner(Corner(triA2, mesh.tris(triA2).c[c]));
+ if (haveB) {
+ for (int c = 0; c < 3; c++)
+ mesh.addCorner(Corner(triB1, mesh.tris(triB1).c[c]));
+ for (int c = 0; c < 3; c++)
+ mesh.addCorner(Corner(triB2, mesh.tris(triB2).c[c]));
+ }
+
+ int baseIdx = 3 * (mesh.numTris() - nt);
+ Corner *cBase = &mesh.corners(baseIdx);
+
+ // set next/prev
+ for (int t = 0; t < nt; t++)
+ for (int c = 0; c < 3; c++) {
+ cBase[t * 3 + c].next = baseIdx + t * 3 + ((c + 1) % 3);
+ cBase[t * 3 + c].prev = baseIdx + t * 3 + ((c + 2) % 3);
+ }
+
+ // set opposites
+ // A1
+ cBase[0].opposite = haveB ? (baseIdx + 9) : -1;
+ cBase[1].opposite = baseIdx + 5;
+ cBase[2].opposite = -1;
+ if (ca_old[2].opposite >= 0) {
+ cBase[2].opposite = ca_old[2].opposite;
+ mesh.corners(cBase[2].opposite).opposite = baseIdx + 2;
+ }
+ // A2
+ cBase[3].opposite = haveB ? (baseIdx + 6) : -1;
+ cBase[4].opposite = -1;
+ if (ca_old[1].opposite >= 0) {
+ cBase[4].opposite = ca_old[1].opposite;
+ mesh.corners(cBase[4].opposite).opposite = baseIdx + 4;
+ }
+ cBase[5].opposite = baseIdx + 1;
+ if (haveB) {
+ // B1
+ cBase[6].opposite = baseIdx + 3;
+ cBase[7].opposite = baseIdx + 11;
+ cBase[8].opposite = -1;
+ if (cb_old[2].opposite >= 0) {
+ cBase[8].opposite = cb_old[2].opposite;
+ mesh.corners(cBase[8].opposite).opposite = baseIdx + 8;
+ }
+ // B2
+ cBase[9].opposite = baseIdx + 0;
+ cBase[10].opposite = -1;
+ if (cb_old[1].opposite >= 0) {
+ cBase[10].opposite = cb_old[1].opposite;
+ mesh.corners(cBase[10].opposite).opposite = baseIdx + 10;
+ }
+ cBase[11].opposite = baseIdx + 7;
+ }
+
+ ////////////////////
+ // mark the two original triangles for deletion
+ taintedTris[triA] = true;
+ mesh.removeTriFromLookup(triA);
+ if (haveB) {
+ taintedTris[triB] = true;
+ mesh.removeTriFromLookup(triB);
+ }
+
+ Real areaA1 = mesh.getFaceArea(triA1), areaA2 = mesh.getFaceArea(triA2);
+ Real areaB1 = 0, areaB2 = 0;
+ if (haveB) {
+ areaB1 = mesh.getFaceArea(triB1);
+ areaB2 = mesh.getFaceArea(triB2);
+ }
+
+ // add channel props for new triangles
+ for (int i = 0; i < mesh.numTriChannels(); i++) {
+ mesh.triChannel(i)->addSplit(triA, areaA1 / (areaA1 + areaA2));
+ mesh.triChannel(i)->addSplit(triA, areaA2 / (areaA1 + areaA2));
+ if (haveB) {
+ mesh.triChannel(i)->addSplit(triB, areaB1 / (areaB1 + areaB2));
+ mesh.triChannel(i)->addSplit(triB, areaB2 / (areaB1 + areaB2));
+ }
+ }
+
+ // add the four new triangles to the prority queue
+ for (int i = mesh.numTris() - nt; i < mesh.numTris(); i++) {
+ // find the maximum length edge in this triangle
+ Vec3 ne0 = mesh.getEdge(i, 0), ne1 = mesh.getEdge(i, 1), ne2 = mesh.getEdge(i, 2);
+ Real nd0 = normSquare(ne0);
+ Real nd1 = normSquare(ne1);
+ Real nd2 = normSquare(ne2);
+ Real longest = max(nd0, max(nd1, nd2));
+ // longest = (int)(longest * 1e2) / 1e2; // HACK: truncate
+ pq.push(pair<Real, int>(longest, i));
+ }
+ edgeSubdivs++;
+ }
+ }
+
+ //////////////////////////////////////////
+ // EDGE COLLAPSING //
+ // - based on short edge length //
+ //////////////////////////////////////////
+ if (minLength > 0) {
+ const Real minLength2 = minLength * minLength;
+ for (int t = 0; t < mesh.numTris(); t++) {
+ // we only want to run through the edge list ONCE.
+ // we achieve this in a method very similar to the above subdivision method.
+
+ // NOTE:
+ // priority queue does not work so great in the edge collapse case,
+ // because collapsing one triangle affects the edge lengths
+ // of many neighbor triangles,
+ // and we do not update their maximum edge length in the queue.
+
+ // if this triangle has already been deleted, ignore it
+ // if(taintedTris[t])
+ // continue;
+
+ if (taintedTris.find(t) != taintedTris.end())
+ continue;
+
+ // first we find the minimum length edge in this triangle
+ Vec3 e0 = mesh.getEdge(t, 0), e1 = mesh.getEdge(t, 1), e2 = mesh.getEdge(t, 2);
+ Real d0 = normSquare(e0);
+ Real d1 = normSquare(e1);
+ Real d2 = normSquare(e2);
+
+ Vec3 edgevect;
+ Vec3 endpoint;
+ Real dist2;
+ int which;
+ if (d0 < d1) {
+ if (d0 < d2) {
+ dist2 = d0;
+ edgevect = e0;
+ endpoint = mesh.getNode(t, 0);
+ which = 2; // 2 opposite of edge 0-1
+ }
+ else {
+ dist2 = d2;
+ edgevect = e2;
+ endpoint = mesh.getNode(t, 2);
+ which = 1; // 1 opposite of edge 2-0
+ }
+ }
+ else {
+ if (d1 < d2) {
+ dist2 = d1;
+ edgevect = e1;
+ endpoint = mesh.getNode(t, 1);
+ which = 0; // 0 opposite of edge 1-2
+ }
+ else {
+ dist2 = d2;
+ edgevect = e2;
+ endpoint = mesh.getNode(t, 2);
+ which = 1; // 1 opposite of edge 2-0
+ }
+ }
+ // then we see if the min length edge is too short
+ if (dist2 < minLength2) {
+ CollapseEdge(
+ mesh, t, which, edgevect, endpoint, deletedNodes, taintedTris, edgeCollsLen, cutTubes);
+ }
+ }
+ }
+ // cleanup nodes and triangles marked for deletion
+
+ // we run backwards through the deleted array,
+ // replacing triangles with ones from the back
+ // (this avoids the potential problem of overwriting a triangle
+ // with a to-be-deleted triangle)
+ std::map<int, bool>::reverse_iterator tti = taintedTris.rbegin();
+ for (; tti != taintedTris.rend(); tti++)
+ mesh.removeTri(tti->first);
+
+ mesh.removeNodes(deletedNodes);
+ cout << "Surface subdivision finished with " << mesh.numNodes() << " surface nodes and "
+ << mesh.numTris();
+ cout << " surface triangles, edgeSubdivs:" << edgeSubdivs << ", edgeCollapses: " << edgeCollsLen;
+ cout << " + " << edgeCollsAngle << " + " << edgeKill << endl;
+ // mesh.sanityCheck();
+}
+static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+{
+ try {
+ PbArgs _args(_linargs, _kwds);
+ FluidSolver *parent = _args.obtainParent();
+ bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+ pbPreparePlugin(parent, "subdivideMesh", !noTiming);
+ PyObject *_retval = 0;
+ {
+ ArgLocker _lock;
+ Mesh &mesh = *_args.getPtr<Mesh>("mesh", 0, &_lock);
+ Real minAngle = _args.get<Real>("minAngle", 1, &_lock);
+ Real minLength = _args.get<Real>("minLength", 2, &_lock);
+ Real maxLength = _args.get<Real>("maxLength", 3, &_lock);
+ bool cutTubes = _args.getOpt<bool>("cutTubes", 4, false, &_lock);
+ _retval = getPyNone();
+ subdivideMesh(mesh, minAngle, minLength, maxLength, cutTubes);
+ _args.check();
+ }
+ pbFinalizePlugin(parent, "subdivideMesh", !noTiming);
+ return _retval;
+ }
+ catch (std::exception &e) {
+ pbSetError("subdivideMesh", e.what());
+ return 0;
+ }
+}
+static const Pb::Register _RP_subdivideMesh("", "subdivideMesh", _W_1);
+extern "C" {
+void PbRegister_subdivideMesh()
+{
+ KEEP_UNUSED(_RP_subdivideMesh);
+}
+}
+
+void killSmallComponents(Mesh &mesh, int elements = 10)
+{
+ const int num = mesh.numTris();
+ vector<int> comp(num);
+ vector<int> numEl;
+ vector<int> deletedNodes;
+ vector<bool> isNodeDel(mesh.numNodes());
+ map<int, bool> taintedTris;
+ // enumerate components
+ int cur = 0;
+ for (int i = 0; i < num; i++) {
+ if (comp[i] == 0) {
+ cur++;
+ comp[i] = cur;
+
+ stack<int> stack;
+ stack.push(i);
+ int cnt = 1;
+ while (!stack.empty()) {
+ int tri = stack.top();
+ stack.pop();
+ for (int c = 0; c < 3; c++) {
+ int op = mesh.corners(tri, c).opposite;
+ if (op < 0)
+ continue;
+ int ntri = mesh.corners(op).tri;
+ if (comp[ntri] == 0) {
+ comp[ntri] = cur;
+ stack.push(ntri);
+ cnt++;
+ }
+ }
+ }
+ numEl.push_back(cnt);
+ }
+ }
+ // kill small components
+ for (int j = 0; j < num; j++) {
+ if (numEl[comp[j] - 1] < elements) {
+ taintedTris[j] = true;
+ for (int c = 0; c < 3; c++) {
+ int n = mesh.tris(j).c[c];
+ if (!isNodeDel[n]) {
+ isNodeDel[n] = true;
+ deletedNodes.push_back(n);
+ }
+ }
+ }
+ }
+
+ std::map<int, bool>::reverse_iterator tti = taintedTris.rbegin();
+ for (; tti != taintedTris.rend(); tti++)
+ mesh.removeTri(tti->first);
+
+ mesh.removeNodes(deletedNodes);
+
+ if (!taintedTris.empty())
+ cout << "Killed small components : " << deletedNodes.size() << " nodes, " << taintedTris.size()
+ << " tris deleted." << endl;
+}
+static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+{
+ try {
+ PbArgs _args(_linargs, _kwds);
+ FluidSolver *parent = _args.obtainParent();
+ bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+ pbPreparePlugin(parent, "killSmallComponents", !noTiming);
+ PyObject *_retval = 0;
+ {
+ ArgLocker _lock;
+ Mesh &mesh = *_args.getPtr<Mesh>("mesh", 0, &_lock);
+ int elements = _args.getOpt<int>("elements", 1, 10, &_lock);
+ _retval = getPyNone();
+ killSmallComponents(mesh, elements);
+ _args.check();
+ }
+ pbFinalizePlugin(parent, "killSmallComponents", !noTiming);
+ return _retval;
+ }
+ catch (std::exception &e) {
+ pbSetError("killSmallComponents", e.what());
+ return 0;
+ }
+}
+static const Pb::Register _RP_killSmallComponents("", "killSmallComponents", _W_2);
+extern "C" {
+void PbRegister_killSmallComponents()
+{
+ KEEP_UNUSED(_RP_killSmallComponents);
+}
+}
+
+} // namespace Manta