------------------------------------------------------------------------------
| 96 | |
| 97 | //------------------------------------------------------------------------------ |
| 98 | int vtkQuad::EvaluatePosition(const double x[3], double closestPoint[3], int& subId, |
| 99 | double pcoords[3], double& dist2, double weights[]) |
| 100 | { |
| 101 | int i, j; |
| 102 | const double *pt1, *pt2, *pt3, *pt; |
| 103 | double n[3]; |
| 104 | double det; |
| 105 | double maxComponent; |
| 106 | int idx = 0, indices[2]; |
| 107 | int iteration, converged; |
| 108 | double params[2]; |
| 109 | double fcol[2], rcol[2], scol[2], cp[3]; |
| 110 | double derivs[8]; |
| 111 | |
| 112 | subId = 0; |
| 113 | pcoords[0] = pcoords[1] = params[0] = params[1] = 0.5; |
| 114 | pcoords[2] = 0.0; |
| 115 | |
| 116 | // Efficient point access |
| 117 | const auto pointsArray = vtkDoubleArray::FastDownCast(this->Points->GetData()); |
| 118 | if (!pointsArray) |
| 119 | { |
| 120 | vtkErrorMacro(<< "Points should be double type"); |
| 121 | return 0; |
| 122 | } |
| 123 | const double* pts = pointsArray->GetPointer(0); |
| 124 | |
| 125 | // Get normal for quadrilateral |
| 126 | // |
| 127 | pt1 = pts; |
| 128 | pt2 = pts + 3; |
| 129 | pt3 = pts + 6; |
| 130 | ComputeNormal(this, pt1, pt2, pt3, n); |
| 131 | |
| 132 | // Project point to plane |
| 133 | // |
| 134 | vtkPlane::ProjectPoint(x, pt1, n, cp); |
| 135 | |
| 136 | // Construct matrices. Since we have over determined system, need to find |
| 137 | // which 2 out of 3 equations to use to develop equations. (Any 2 should |
| 138 | // work since we've projected point to plane.) |
| 139 | // |
| 140 | for (maxComponent = 0.0, i = 0; i < 3; i++) |
| 141 | { |
| 142 | if (fabs(n[i]) > maxComponent) |
| 143 | { |
| 144 | maxComponent = fabs(n[i]); |
| 145 | idx = i; |
| 146 | } |
| 147 | } |
| 148 | for (j = 0, i = 0; i < 3; i++) |
| 149 | { |
| 150 | if (i != idx) |
| 151 | { |
| 152 | indices[j++] = i; |
| 153 | } |
| 154 | } |
| 155 |