| 94 | constexpr float eps = 0.001f; |
| 95 | |
| 96 | void updateIndicator( const Mesh & mesh, Vector<float, UndirectedEdgeId> & v, const FaceNormals & normals, float beta, float gamma ) |
| 97 | { |
| 98 | MR_TIMER; |
| 99 | |
| 100 | const auto sz = v.size(); |
| 101 | assert( sz == mesh.topology.undirectedEdgeSize() ); |
| 102 | assert( (int)normals.size() >= mesh.topology.lastValidFace() ); |
| 103 | if ( sz <= 0 ) |
| 104 | return; |
| 105 | |
| 106 | std::vector< Eigen::Triplet<double> > mTriplets; |
| 107 | Eigen::VectorXd rhs; |
| 108 | rhs.resize( sz ); |
| 109 | const float rh = beta / ( 2 * eps ); |
| 110 | const float k = 2 * beta * eps; |
| 111 | for ( auto ue = 0_ue; ue < sz; ++ue ) |
| 112 | { |
| 113 | const EdgeId e = ue; // note that it can be lone edge |
| 114 | float centralWeight = rh; |
| 115 | const auto l = mesh.topology.left( e ); |
| 116 | const auto r = mesh.topology.right( e ); |
| 117 | if ( l && r ) |
| 118 | centralWeight += 2 * gamma * ( normals[l] - normals[r] ).lengthSq(); |
| 119 | const auto lenE = ( l || r ) ? mesh.edgeLength( e ) : 0.0f; |
| 120 | if ( lenE > 0 ) |
| 121 | { |
| 122 | if ( l ) |
| 123 | { |
| 124 | const auto c = mesh.triCenter( l ); |
| 125 | { |
| 126 | const auto a = mesh.topology.next( e ); |
| 127 | const auto lenL = ( c - mesh.orgPnt( e ) ).length(); |
| 128 | const auto x = k * lenL / lenE; |
| 129 | centralWeight += x; |
| 130 | mTriplets.emplace_back( ue, a.undirected(), -x ); |
| 131 | } |
| 132 | { |
| 133 | const auto b = mesh.topology.prev( e.sym() ); |
| 134 | const auto lenL = ( c - mesh.destPnt( e ) ).length(); |
| 135 | const auto x = k * lenL / lenE; |
| 136 | centralWeight += x; |
| 137 | mTriplets.emplace_back( ue, b.undirected(), -x ); |
| 138 | } |
| 139 | } |
| 140 | if ( r ) |
| 141 | { |
| 142 | const auto c = mesh.triCenter( r ); |
| 143 | { |
| 144 | const auto a = mesh.topology.prev( e ); |
| 145 | const auto lenL = ( c - mesh.orgPnt( e ) ).length(); |
| 146 | const auto x = k * lenL / lenE; |
| 147 | centralWeight += x; |
| 148 | mTriplets.emplace_back( ue, a.undirected(), -x ); |
| 149 | } |
| 150 | { |
| 151 | const auto b = mesh.topology.next( e.sym() ); |
| 152 | const auto lenL = ( c - mesh.destPnt( e ) ).length(); |
| 153 | const auto x = k * lenL / lenE; |
no test coverage detected