Skip to content
Snippets Groups Projects
Commit b4c92cf0 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

*** empty log message ***

parent 8797da85
No related branches found
No related tags found
No related merge requests found
......@@ -190,53 +190,35 @@ GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2])
double XP = qp.x();
double YP = qp.y();
double ZP = qp.z();
//double a = meanPlane.a;
//double b = meanPlane.b;
//double c = meanPlane.c;
//double d = meanPlane.d;
// X = P + H N, find y such that
// a X_1 + b X_2 + c X_3 + d = 0
// a ( XP + y XN) + b ( YP + y YN) + c ( ZP + y ZN) + d = 0
// H ( a XN + b YN + c ZN ) = - (a XP + b YP + c ZP + d)
//const double H = -(a*XP + b*YP + c *ZP + d)/(a*a + b*b + c*c);
//const double X = XP + H * a;
//const double Y = YP + H * b;
//const double Z = ZP + H * c;
// now compute parametric coordinates
double x,y,z,VX[3], VY[3];
double VX[3], VY[3], x, y, z;
getMeanPlaneData(VX, VY, x, y, z);
// XP = X + u VX + v VY
// We are sure to be on the plane
double M[3][2] = {{VX[0],VY[0]},{VX[1],VY[1]},{VX[2],VY[2]}};
double M[3][2] = {{VX[0], VY[0]}, {VX[1], VY[1]}, {VX[2], VY[2]}};
double MN[2][2];
double B[3] = {XP-x,YP-y,ZP-z};
double BN[2],UV[2];
for (int i=0;i<2;i++){
double B[3] = {XP - x, YP - y, ZP - z};
double BN[2], UV[2];
for (int i = 0; i < 2; i++){
BN[i] = 0;
for (int k=0;k<3;k++){
for (int k = 0; k < 3; k++){
BN[i] += B[k] * M[k][i];
}
}
for (int i=0;i<2;i++){
for (int j=0;j<2;j++){
for (int i = 0; i < 2; i++){
for (int j = 0; j < 2; j++){
MN[i][j] = 0;
for (int k=0;k<3;k++){
for (int k = 0; k < 3; k++){
MN[i][j] += M[k][i] * M[k][j];
}
}
}
sys2x2(MN,BN,UV);
// GPoint test = point (UV[0],UV[1]);
// printf("%g %g %g vs %g %g %g\n",XP,YP,ZP,test.x(),test.y(),test.z());
sys2x2(MN, BN, UV);
return GPoint(XP, YP, ZP, this, UV);
}
Vertex v;
double u[2] = {initialGuess[0],initialGuess[1]};
v.Pos.X = qp.x();
v.Pos.Y = qp.y();
v.Pos.Z = qp.z();
double u[2] = {initialGuess[0], initialGuess[1]};
bool result = ProjectPointOnSurface(s, v, u);
if (!result)
return GPoint(-1.e22, -1.e22, -1.e22, 0, u);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment