| 92 | } |
| 93 | |
| 94 | void TetModel::attachVisMesh(const ParticleData &pd) |
| 95 | { |
| 96 | const Real eps = static_cast<Real>(1.0e-6); |
| 97 | |
| 98 | // The created surface mesh defines the boundary of the tet mesh |
| 99 | unsigned int *faces = m_surfaceMesh.getFaces().data(); |
| 100 | const unsigned int nFaces = m_surfaceMesh.numFaces(); |
| 101 | |
| 102 | const Vector3r *normals = m_surfaceMesh.getVertexNormals().data(); |
| 103 | |
| 104 | // for each point find nearest triangle (TODO: optimize) |
| 105 | const int nNearstT = 15; |
| 106 | m_attachments.resize(m_visVertices.size()); |
| 107 | |
| 108 | #pragma omp parallel default(shared) |
| 109 | { |
| 110 | #pragma omp for schedule(static) |
| 111 | for (int i = 0; i < (int)m_visVertices.size(); i++) |
| 112 | { |
| 113 | const Vector3r &p = m_visVertices.getPosition(i); |
| 114 | Real curDist[nNearstT]; |
| 115 | int curT[nNearstT]; |
| 116 | for (int k = 0; k < nNearstT; k++) |
| 117 | { |
| 118 | curDist[k] = REAL_MAX; |
| 119 | curT[k] = -1; |
| 120 | } |
| 121 | Vector3r curBary[nNearstT]; |
| 122 | Vector3r curInter[nNearstT]; |
| 123 | for (unsigned int j = 0; j < nFaces; j++) |
| 124 | { |
| 125 | const unsigned int indexA = faces[3 * j] + m_indexOffset; |
| 126 | const unsigned int indexB = faces[3 * j + 1] + m_indexOffset; |
| 127 | const unsigned int indexC = faces[3 * j + 2] + m_indexOffset; |
| 128 | const Vector3r & a = pd.getPosition0(indexA); |
| 129 | const Vector3r & b = pd.getPosition0(indexB); |
| 130 | const Vector3r & c = pd.getPosition0(indexC); |
| 131 | |
| 132 | Vector3r inter, bary; |
| 133 | // compute nearest point on triangle |
| 134 | if (pointInTriangle(a, b, c, p, inter, bary)) |
| 135 | { |
| 136 | Real len = (p - inter).norm(); |
| 137 | for (int k = nNearstT - 1; k >= 0; k--) // update the best triangles |
| 138 | { |
| 139 | if (len < curDist[k]) |
| 140 | { |
| 141 | if (k < nNearstT - 1) |
| 142 | { |
| 143 | curDist[k + 1] = curDist[k]; |
| 144 | curBary[k + 1] = curBary[k]; |
| 145 | curT[k + 1] = curT[k]; |
| 146 | curInter[k + 1] = curInter[k]; |
| 147 | } |
| 148 | curDist[k] = len; |
| 149 | curBary[k] = bary; |
| 150 | curT[k] = (int)j; |
| 151 | curInter[k] = inter; |